diff --git a/Data/proteinAllocationModel_iML1515_EnzymaticData_py_uniprot.xlsx b/Data/proteinAllocationModel_iML1515_EnzymaticData_py_uniprot.xlsx
new file mode 100644
index 0000000..be99597
Binary files /dev/null and b/Data/proteinAllocationModel_iML1515_EnzymaticData_py_uniprot.xlsx differ
diff --git a/Results/protein_costs_glycolysis_tca.xlsx b/Results/protein_costs_glycolysis_tca.xlsx
new file mode 100644
index 0000000..a9cd99a
Binary files /dev/null and b/Results/protein_costs_glycolysis_tca.xlsx differ
diff --git a/Scripts/pam_generation_new.py b/Scripts/pam_generation_new.py
new file mode 100644
index 0000000..989bc0b
--- /dev/null
+++ b/Scripts/pam_generation_new.py
@@ -0,0 +1,452 @@
+import cobra
+import pandas as pd
+import numpy as np
+import os
+from typing import Union
+
+# load PAMpy modules
+from src.PAModelpy.PAModel import PAModel
+from src.PAModelpy.EnzymeSectors import ActiveEnzymeSector, UnusedEnzymeSector, TransEnzymeSector
+from src.PAModelpy.configuration import Config
+
+from Scripts.toy_ec_pam import build_toy_gem, build_active_enzyme_sector, build_translational_protein_sector, build_unused_protein_sector
+import ast
+
+'Function library for making Protein Allocation Models as described in the publication'
+
+
+def set_up_toy_pam(sensitivity =True):
+    config = Config()
+    #setting the configuration for the toy model
+    config.BIOMASS_REACTION = 'R7'
+    config.GLUCOSE_EXCHANGE_RXNID = 'R1'
+    config.CO2_EXHANGE_RXNID = 'R8'
+    config.ACETATE_EXCRETION_RXNID = 'R9'
+
+    Etot = 0.6*1e-3
+    model = build_toy_gem()
+    active_enzyme = build_active_enzyme_sector(config)
+    unused_enzyme = build_unused_protein_sector(config)
+    translation_enzyme = build_translational_protein_sector(config)
+    pamodel = PAModel(model, name='toy model MCA with enzyme constraints', active_sector=active_enzyme,
+                      translational_sector=translation_enzyme,
+                      unused_sector=unused_enzyme, p_tot=Etot, sensitivity=sensitivity)
+    pamodel.objective = 'R7'
+    config.reset()
+    return pamodel
+
+def set_up_ecolicore_pam(total_protein:bool = True, active_enzymes: bool = True, translational_enzymes:bool = True, unused_enzymes:bool = True, sensitivity =True):
+    # Setting the relative paths
+    PAM_DATA_FILE_PATH = os.path.join('Data', 'proteinAllocationModel_iML1515_EnzymaticData_py_uniprot.xlsx')
+
+    config = Config()
+    config.reset()
+
+    # some other constants
+    BIOMASS_REACTION = 'BIOMASS_Ecoli_core_w_GAM'
+    TOTAL_PROTEIN_CONCENTRATION = 0.16995  # [g_prot/g_cdw]
+
+    # load the genome-scale information
+    model = cobra.io.load_json_model(os.path.join('Models', 'e_coli_core.json'))
+
+    #load example data for the E.coli iML1515 model
+    if active_enzymes:
+        # load active enzyme sector information
+        enzyme_db = pd.read_excel(PAM_DATA_FILE_PATH, sheet_name='ActiveEnzymes')
+        # create enzyme objects for each gene-associated reaction
+        rxn2protein, protein2gene = parse_reaction2protein(enzyme_db, model)
+
+        rxn2protein_new = []
+        for i in list(rxn2protein.keys())[:5]:
+            print(i)
+            rxn2protein_new[i] = rxn2protein[i].copy()
+
+        rxn2protein = rxn2protein_new
+        active_enzyme_sector = ActiveEnzymeSector(rxn2protein=rxn2protein, protein2gene=protein2gene,
+                                                  configuration=config)
+
+    else:
+        active_enzyme_sector = None
+
+    if translational_enzymes:
+        # translational protein sector parameter (substrate dependent)
+        id_list_tps = ['EX_glc__D_e']
+        tps_0 = [0.04992]  # g/gDW
+        tps_mu = [-0.002944]  # g h/gDW -> transformed to match glucose uptake variable
+        molmass_tps = [405903.94]  # g/mol
+
+        # translational protein sector
+        translation_enzyme_sector = TransEnzymeSector(
+            id_list=id_list_tps,
+            tps_0=tps_0,
+            tps_mu=tps_mu,
+            mol_mass=molmass_tps,
+        )
+    else:
+        translation_enzyme_sector = None
+
+    if unused_enzymes:
+        id_list_ups = [BIOMASS_REACTION]
+        ups_0 = [0.0407]  # g/gDW
+        ups_mu = [-0.0214]  # g h/gDW -> negative relation with growth rate
+        molmass_ups = [405903.94]  # g/mol
+
+        unused_enzyme_sector = UnusedEnzymeSector(
+            id_list=id_list_ups,
+            ups_0=ups_0,
+            ups_mu=ups_mu,
+            mol_mass=molmass_ups,
+        )
+    else:
+        unused_enzyme_sector = None
+
+    if total_protein: total_protein = TOTAL_PROTEIN_CONCENTRATION
+
+    pa_model = PAModel(id_or_model=model, p_tot=total_protein, sensitivity=sensitivity,
+                       active_sector=active_enzyme_sector, translational_sector=translation_enzyme_sector, unused_sector=unused_enzyme_sector)
+    return pa_model
+
+def set_up_ecoli_pam(total_protein: Union[bool, float] = True, active_enzymes: bool = True,
+                   translational_enzymes: bool = True, unused_enzymes: bool = True, sensitivity = True):
+
+    config = Config()
+    config.reset()
+
+    pam_info_file = os.path.join('Data', 'proteinAllocationModel_iML1515_EnzymaticData_py_uniprot.xlsx')
+
+    # some other constants
+    TOTAL_PROTEIN_CONCENTRATION = 0.258  # [g_prot/g_cdw]
+
+    #setup the gem ecoli iML1515 model
+    model = cobra.io.read_sbml_model(os.path.join('Models', 'iML1515.xml'))
+
+    #check if a different total protein concentration is given
+    if isinstance(total_protein, float):
+        TOTAL_PROTEIN_CONCENTRATION = total_protein
+
+    # load example data for the E.coli iML1515 model
+    if active_enzymes:
+        # load active enzyme sector information
+        enzyme_db = pd.read_excel(pam_info_file, sheet_name='ActiveEnzymes')
+
+        # create enzyme objects for each gene-associated reaction
+        rxn2protein, protein2gene = parse_reaction2protein(enzyme_db, model)
+
+        active_enzyme_sector = ActiveEnzymeSector(rxn2protein=rxn2protein, protein2gene=protein2gene,
+                                                  configuration=config)
+
+    else:
+        active_enzyme_sector = None
+
+    if translational_enzymes:
+        translational_info = pd.read_excel(pam_info_file, sheet_name='Translational')
+        translation_enzyme_sector = TransEnzymeSector(
+            id_list=[translational_info[translational_info.Parameter == 'id_list'].loc[0, 'Value']],
+            tps_0=[translational_info[translational_info.Parameter == 'tps_0'].loc[1, 'Value']],
+            tps_mu=[-translational_info[translational_info.Parameter == 'tps_mu'].loc[2, 'Value']],
+            mol_mass=[translational_info[translational_info.Parameter == 'mol_mass'].loc[3, 'Value']],
+            configuration = config)
+    else:
+        translation_enzyme_sector = None
+
+    if unused_enzymes:
+        unused_protein_info = pd.read_excel(pam_info_file, sheet_name='ExcessEnzymes')
+
+        ups_0 = unused_protein_info[unused_protein_info.Parameter == 'ups_0'].loc[2, 'Value']
+        smax = unused_protein_info[unused_protein_info.Parameter == 's_max_uptake'].loc[1, 'Value']
+
+        unused_protein_sector = UnusedEnzymeSector(
+            id_list=[unused_protein_info[unused_protein_info.Parameter == 'id_list'].loc[0, 'Value']],
+            ups_mu=[ups_0 / smax],
+            ups_0=[ups_0],
+            mol_mass=[unused_protein_info[unused_protein_info.Parameter == 'mol_mass'].loc[3, 'Value']],
+            configuration = config)
+    else:
+        unused_protein_sector = None
+
+    if total_protein: total_protein = TOTAL_PROTEIN_CONCENTRATION
+
+    pamodel = PAModel(id_or_model=model, p_tot=total_protein,
+                       active_sector=active_enzyme_sector, translational_sector=translation_enzyme_sector,
+                       unused_sector=unused_protein_sector, sensitivity=sensitivity, configuration = config
+                      )
+    return pamodel
+
+
+def _parse_enzyme_information_from_file(file_path:str):
+    # load active enzyme sector information
+    enzyme_db = pd.read_excel(file_path, sheet_name='enzyme-gene-reaction')
+    # enzyme_db = enzyme_db.set_index('rxnID')
+    # correct reaction IDS
+    for idx in enzyme_db.rxnID.to_list():
+        # transport reactions
+        if 'pp' in idx:
+            idx_new = idx.replace('pp', '')
+            if idx_new not in enzyme_db.index:
+                enzyme_db.rename(index={idx: idx_new}, inplace=True)
+        if 'ex' in idx:
+            idx_new = idx.replace('ex', '')
+            if idx_new not in enzyme_db.index:
+                enzyme_db.rename(index={idx: idx_new}, inplace=True)
+
+                # replace NaN values with unique identifiers
+    # replace NaN enzyme ids with a dummy enzyme identifier
+    # select the NaN values
+    nan_values = enzyme_db['EC_nmbr'].isnull()
+    # make a list with unique ids
+    nan_ids = [f'E{i}' for i in range(nan_values.sum())]
+    # replace nan values by unique id
+    enzyme_db.loc[nan_values, 'EC_nmbr'] = nan_ids
+
+    return enzyme_db
+
+
+def parse_reaction2protein(enzyme_db: pd.DataFrame, model:cobra.Model) -> dict:
+    # Initialize dictionaries
+    rxn2protein = {}
+    protein2gpr = {}
+
+    # replace NaN values with unique identifiers
+    # select the NaN values
+    nan_values = enzyme_db['uniprotID'].isnull()
+    # make a list with unique ids
+    nan_ids = [f'E{i}' for i in range(nan_values.sum())]
+    # replace nan values by unique id
+    enzyme_db.loc[nan_values, 'uniprotID'] = nan_ids
+
+    protein2gene, gene2protein = _get_genes_for_proteins(enzyme_db, model)
+
+    # Iterate over each row in the DataFrame
+    for index, row in enzyme_db.iterrows():
+        # Parse data from the row
+        rxn_id = row['rxnID']
+        if rxn_id not in model.reactions: continue
+        # only parse those reactions which are in the model
+        kcat_f_b = [row['kcat_f'], row['kcat_b']]
+        kcat_f_b = [kcat if not np.isnan(kcat) else 0 for kcat in kcat_f_b]
+        if all([np.isnan(kcat) for kcat in kcat_f_b]): continue
+
+        rxn = model.reactions.get_by_id(rxn_id)
+        # get the identifiers and replace nan values by dummy placeholders
+        enzyme_id = row['uniprotID']
+        gene_id = row['m_gene']
+
+        # check if there are genes associates with the reaction
+        if len(rxn.genes) > 0 or isinstance(gene_id, str):
+            if not isinstance(enzyme_id, str):
+                enzyme_id = 'Enzyme_' + rxn_id
+                row['molMass'] = 39959.4825  # default molmass
+            if not isinstance(gene_id, str):
+                gene_id = [gene.id for gene in rxn.genes][0]  # TODO
+
+            # get the gene-protein-reaction-associations for this specific enzyme
+            gr, pr = parse_gpr_information_for_rxn2protein(row['m_gene_reaction_rule'],
+                                                           gene2protein, protein2gene,enzyme_id)
+            if pr is None: pr = [[enzyme_id]]
+
+            if enzyme_id not in protein2gpr:
+                protein2gpr[enzyme_id] = gr
+            else:
+                protein2gpr[enzyme_id].append(gr)
+
+            # Create rxn2protein dictionary
+            if rxn_id not in rxn2protein:
+                rxn2protein[rxn_id] = {}
+            if enzyme_id not in rxn2protein[rxn_id]:
+                rxn2protein[rxn_id][enzyme_id] = {
+                    'f': kcat_f_b[0],  # Forward kcat
+                    'b': kcat_f_b[1],  # Backward kcat
+                    'molmass': row['molMass'],
+                    'genes': gr,
+                    'protein_reaction_association': pr,
+                    'complex_with': None
+                }
+            else:
+                rxn2protein[rxn_id][enzyme_id]['genes'].append(gene_id)
+
+    # if no enzyme info is found, add dummy enzyme with median kcat and molmass
+    for rxn in model.reactions:
+        if rxn.id not in rxn2protein.keys() and 'EX'.lower() not in rxn.id.lower() and 'BIOMASS' not in rxn.id and len(
+                rxn._genes) > 0 and list(rxn._genes)[0].id != 's0001':
+            rev = _check_reaction_reversibility(rxn)
+            if rev == 0:
+                kcat_dict = {'f': 22}
+            elif rev == 1:
+                kcat_dict = {'b': 22}
+            else:
+                kcat_dict = {'f': 22, 'b': 22}
+            # no enzyme information found
+            print('No enzyme information found for reaction: ' + rxn.id)
+            enzyme_id = 'Enzyme_' + rxn.id
+            gpr_info = parse_gpr_information_for_protein2genes(rxn.gpr, rxn_id)
+
+            rxn2protein[rxn.id] = {enzyme_id: {
+                **kcat_dict,
+                'molmass': 3.947778784340140e04,
+                'genes': gpr_info
+            }}
+            # add geneinfo for unknown enzymes
+            protein2gpr[enzyme_id] = gpr_info
+
+        return rxn2protein, protein2gpr
+
+def _get_genes_for_proteins(enzyme_db: pd.DataFrame, model) -> dict:
+    protein2gene = {}
+    gene2protein = {}
+    for index, row in enzyme_db.iterrows():
+        # Parse data from the row
+        rxn_id = row['rxnID']
+        if rxn_id not in model.reactions:continue
+        rxn = model.reactions.get_by_id(rxn_id)
+        # get the identifiers and replace nan values by dummy placeholders
+        enzyme_id = row['uniprotID']
+        gene_id = row['m_gene']
+
+        # check if there are genes associates with the reaction
+        if len(rxn.genes) > 0 or isinstance(gene_id, str):
+            if not isinstance(enzyme_id, str):
+                enzyme_id = 'Enzyme_' + rxn_id
+                row['molMass'] = 39959.4825  # default molmass
+            if not isinstance(gene_id, str):
+                gene_id = [gene.id for gene in rxn.genes][0]  # TODO
+
+            gene2protein[gene_id] = enzyme_id
+
+            # Create gene-protein-reaction associations
+            if enzyme_id not in protein2gene:
+                protein2gene[enzyme_id] = [gene_id]
+            elif enzyme_id in protein2gene:
+                # assume that there is a single protein if the previous relation was not assigned a gene
+                if 'gene_' in protein2gene[enzyme_id][0]:
+                    protein2gene[enzyme_id] = [gene_id]
+                else:
+                    protein2gene[enzyme_id].append(gene_id)
+
+    return protein2gene, gene2protein
+
+
+def _check_rxn_identifier_format(rxn_id:str) -> str:
+    if 'pp' in rxn_id:
+        idx_new = rxn_id.replace('pp', '')
+    elif 'ex' in rxn_id:
+        idx_new = rxn_id.replace('ex', '')
+    else:
+        idx_new = rxn_id
+    return idx_new
+
+def _get_fwd_bckw_kcat(rxn_id: str, kcat:float, model:PAModel) -> Union[list, None]:
+    #skip exchange and biomass reaction
+    if 'EX' in rxn_id or 'BIOMASS' in rxn_id:
+        return None
+
+    # Extract the base identifier without any suffixes
+    base_id = rxn_id.split('_')[0]
+
+    # Iterate over each identifier in the input
+    if base_id in model.reactions:
+        # Determine the form of the identifier
+        if rxn_id.endswith('_f'):
+            kcat_fwd = kcat
+            kcat_rev = 0
+        elif rxn_id.endswith('_b'):
+            kcat_fwd = 0
+            kcat_rev = kcat
+        # identifier has no suffix
+        elif base_id == rxn_id:
+            kcat_fwd = kcat
+            kcat_rev = kcat
+        else:
+            return None
+    elif rxn_id in model.reactions:
+        kcat_fwd = kcat
+        kcat_rev = kcat
+    else:
+        return None
+    return [kcat_fwd, kcat_rev]
+
+def _check_reaction_reversibility(reaction):
+    if reaction.lower_bound >= 0:
+        # irreversible reaction (forward direction)
+        rev = 0
+    elif reaction.upper_bound <= 0:
+        # irreversible reaction (reverse direction)
+        rev = 1
+    else:
+        rev = 2
+        # reversible r
+    return rev
+
+def parse_gpr_information_for_protein2genes(gpr_info:str):
+    #filter out nan entries
+    if isinstance(gpr_info, str):
+        gpr_list = parse_gpr(gpr_info)
+        # Extracting the inner lists and removing parentheses
+        # gpr_list = [[[item.strip("(')")] + sublist[1:] for item in sublist[0].split(" or ")] for sublist in nested_list]
+        return gpr_list
+    else:
+        return [['gene_dummy']]
+
+def parse_gpr_information_for_rxn2protein(gpr_info:str, gene2protein: dict,protein2gene:dict, enzyme_id:str):
+    #filter out nan entries
+    if isinstance(gpr_info, str):
+        gpr_list = parse_gpr(gpr_info)
+        # Extracting the inner lists and removing parentheses
+        gene = protein2gene[enzyme_id]
+        if isinstance(protein2gene[enzyme_id], list):
+            gene = gene[0]
+        gpr_list = filter_sublists(gpr_list, gene)
+        #convert the genes to the associated proteins
+        enzyme_relations = []
+        for sublist in gpr_list:
+            enz_sublist = []
+            for item in sublist:
+                if item in gene2protein.keys():
+                    enz_sublist.append(gene2protein[item])
+            enzyme_relations += [enz_sublist]
+        enzyme_relations = filter_sublists(enzyme_relations, enzyme_id)
+        return gpr_list, enzyme_relations
+    else:
+        return [['gene_dummy']], None
+
+def parse_gpr(gpr_info):
+    # Split the string by 'or' first
+    or_groups = gpr_info.split(' or ')
+    parsed_gpr = []
+
+    for group in or_groups:
+        # Within each 'or' group, split by 'and'
+        and_genes = group.split(' and ')
+        # Clean up each gene string
+        and_genes = [gene.strip(' ()') for gene in and_genes]
+        parsed_gpr.append(and_genes)
+
+    return parsed_gpr
+
+def filter_sublists(nested_list, target_string):
+    """
+    Filters out all sublists from a nested list that contain the target string.
+
+    Args:
+        nested_list (list of list of str): The nested list to filter.
+        target_string (str): The string to filter out sublists that contain it.
+
+    Returns:
+        list of list of str: A new nested list with the filtered sublists.
+    """
+    return [sublist for sublist in nested_list if target_string in sublist]
+
+if __name__ == '__main__':
+    ecoli_pam = set_up_ecolicore_pam(sensitivity=False)
+    # ecoli_pam.objective = ecoli_pam.BIOMASS_REACTION
+    ecoli_pam.change_reaction_bounds('EX_glc__D_e', -10, 0)
+    ecoli_pam.optimize()
+    print(ecoli_pam.objective.value)
+
+    for enzyme in ecoli_pam.enzymes:
+        if hasattr(enzyme, 'enzymes'):
+            print(enzyme)
+            for c in enzyme._constraints.values():
+                print(c)
+    for ce in ecoli_pam.catalytic_events.query('AKGDH'):
+        for c in ce.constraints.values():
+            print(c)
diff --git a/Scripts/protein_costs_analysis.py b/Scripts/protein_costs_analysis.py
index 7db7e06..63a1da8 100644
--- a/Scripts/protein_costs_analysis.py
+++ b/Scripts/protein_costs_analysis.py
@@ -1,9 +1,8 @@
 import os
 import pandas as pd
 
-DATA_DIR = os.path.join(os.path.split(os.getcwd())[0], 'Data')
-PAM_DATA_FILE_PATH = os.path.join(DATA_DIR, 'proteinAllocationModel_iML1515_EnzymaticData_py.xls')
-RESULT_DATA_FILE = os.path.join(os.path.split(os.getcwd())[0], 'Results', 'protein_costs_glycolysis_tca.xlsx')
+PAM_DATA_FILE_PATH = os.path.join('Data', 'proteinAllocationModel_iML1515_EnzymaticData_py.xls')
+RESULT_DATA_FILE = os.path.join('Results', 'protein_costs_glycolysis_tca.xlsx')
 
 active_enzyme_info = pd.read_excel(PAM_DATA_FILE_PATH, sheet_name='ActiveEnzymes')
 
@@ -15,7 +14,9 @@
 protein_efficiency_df = active_enzyme_info.filter(['rxnID','molMass','kcat'])
 protein_efficiency_df = protein_efficiency_df.assign(rxnName = lambda x: x.rxnID.str.rstrip('_fb'))
 protein_efficiency_glyc_tca = protein_efficiency_df[protein_efficiency_df.rxnName.isin(tca_glycolysis_reactions)]
-protein_efficiency_glyc_tca  = protein_efficiency_glyc_tca.assign(relEfficiency = lambda x: x['molMass']/x['kcat']).sort_values(by = 'relEfficiency', ascending=False)
+protein_efficiency_glyc_tca  = protein_efficiency_glyc_tca.assign(kcat = lambda x:x['kcat']*3600,
+                                                                  relEfficiency = lambda x: x['molMass']/x['kcat']
+                                                                  ).sort_values(by = 'relEfficiency', ascending=False)
 print(protein_efficiency_glyc_tca)
 
 protein_efficiency_glyc_tca.to_excel(RESULT_DATA_FILE, sheet_name='protein_costs')
\ No newline at end of file
diff --git a/Scripts/toy_ec_pam.py b/Scripts/toy_ec_pam.py
index de7fe0e..d4c36b8 100644
--- a/Scripts/toy_ec_pam.py
+++ b/Scripts/toy_ec_pam.py
@@ -4,9 +4,9 @@
 import numpy as np
 
 #importing the tools from the PAModelpy package
-from PAModelpy.EnzymeSectors import ActiveEnzymeSector, UnusedEnzymeSector, TransEnzymeSector
-from PAModelpy.PAModel import PAModel
-from PAModelpy.configuration import Config
+from src.PAModelpy.EnzymeSectors import ActiveEnzymeSector, UnusedEnzymeSector, TransEnzymeSector
+from src.PAModelpy.PAModel import PAModel
+from src.PAModelpy.configuration import Config
 
 Config.BIOMASS_REACTION = 'R7'
 Config.GLUCOSE_EXCHANGE_RXNID = 'R1'
@@ -104,14 +104,21 @@ def build_toy_gem():
 def build_active_enzyme_sector(Config):
     kcat_fwd = [1, 0.5, 1, 1, 0.5 ,0.45, 1.5]  # High turnover for exchange reactions
     kcat_rev = [kcat for kcat in kcat_fwd]
+    protein2gene = {}
     rxn2kcat = {}
     for i in range(n-3): # all reactions have an enzyme, except excretion reactions
         rxn_id = f'R{i+1}'
         # 1e-6 to correct for the unit transformation in the model (meant to make the calculations preciser for different unit dimensions)
         #dummy molmass like in MATLAB script
-        rxn2kcat = {**rxn2kcat, **{rxn_id: {f'E{i+1}':{'f': kcat_fwd[i]/(3600*1e-6), 'b': kcat_rev[i]/(3600*1e-6), 'molmass': 1e6}}}}
-
-    return ActiveEnzymeSector(rxn2protein = rxn2kcat, configuration=Config)
+        rxn2kcat = {**rxn2kcat, **{rxn_id:
+                                       {f'E{i+1}':
+                                            {'f': kcat_fwd[i]/(3600*1e-6),
+                                             'b': kcat_rev[i]/(3600*1e-6),
+                                             'molmass': 1e6,
+                                             'protein_reaction_association':[[f'E{i+1}']]}}}}
+        protein2gene = {**protein2gene,**{f'E{i+1}':[[f'g{i+1}']]}}
+
+    return ActiveEnzymeSector(rxn2protein = rxn2kcat, protein2gene = protein2gene, configuration=Config)
 
 def build_unused_protein_sector(Config):
     return UnusedEnzymeSector(id_list = ['R1'], ups_mu=[-0.01*1e-3], ups_0=[0.1*1e-3], mol_mass= [1], configuration=Config)
diff --git a/src/PAModelpy/CatalyticEvent.py b/src/PAModelpy/CatalyticEvent.py
index 29bf765..c9fbdb9 100644
--- a/src/PAModelpy/CatalyticEvent.py
+++ b/src/PAModelpy/CatalyticEvent.py
@@ -4,14 +4,15 @@
 (e.g. one conversion of substrate to product, can be catalyzed by multiple enzymes)
 """
 import cobra
-from cobra import DictList, Object
+from cobra import DictList, Object, Reaction
 from cobra.exceptions import OptimizationError
 from cobra.util.solver import check_solver_status
 from optlang.symbolics import Zero
-from typing import Optional, Dict
+from typing import Optional, Dict, Union
 from warnings import warn
 from copy import copy, deepcopy
 
+
 class CatalyticEvent(Object):
     """
         CatalyticEvent is a class for holding information regarding the
@@ -49,6 +50,7 @@ def __init__(
         #relation to reaction
         self.rxn_id = rxn_id
         self.rxn = None
+        self.catalytic_reactions = DictList()
 
         #relation to enzymes
         self.kcats = kcats2enzymes
@@ -155,11 +157,22 @@ def model(self, model):
             rxn = cobra.Reaction(id = self.rxn_id)
             self._model.add_reactions([rxn])
             self.rxn = rxn
+        #add reaction constraint
+        ce_constraint = self._model.problem.Constraint(Zero, name= self.id, lb = 0, ub=0)
+        self._model.add_cons_vars([ce_constraint])
+        self._model.constraints[self.id].set_linear_coefficients({
+            self.rxn.forward_variable: 1,
+            self.rxn.reverse_variable: -1,
+        })
+        self.constraints[self.id] = ce_constraint
+
         #add enzymes to the model if they are not there already
         for enzyme in self.enzymes:
             if enzyme in self._model.enzymes:
                 self.constraints ={**self.constraints, **enzyme._constraints}
+                enzyme._constraints[self.id] = ce_constraint
                 enzyme_model = self._model.enzymes.get_by_id(enzyme.id)
+
                 if self.rxn_id in enzyme_model.rxn2kcat.keys():
                     if self not in enzyme_model.catalytic_events:
                         enzyme_model.add_catalytic_event(self, kcats= {enzyme: enzyme.rxn2kcat[self.rxn_id]})
@@ -168,6 +181,11 @@ def model(self, model):
             else:
                 self._model.add_enzymes([enzyme])
 
+            # add enzyme-reaction association
+            self.add_enzyme_reaction_association(enzyme)
+
+
+
     def add_enzymes(self,enzyme_kcat_dict: dict):
         """
         Add enzymes to the catalytic event and create bindings to the related model.
@@ -184,10 +202,8 @@ def add_enzymes(self,enzyme_kcat_dict: dict):
         for enzyme, kcat in enzyme_kcat_dict.items():
             # check if the enzyme is already associated to the catalytic event
             if enzyme in self.enzymes:
-                # print(enzyme)
-                # warn(f'Enzyme {enzyme.id} is already associated with catalytic event {self.id}. This enzyme will be updated')
-                self.change_kcat_values({enzyme.id: enzyme.get_kcat_values(rxn_ids = [self.rxn_id])})
-                # print(self.enzymes)
+                self.add_enzyme_reaction_association(enzyme)
+                self.change_kcat_values({enzyme.id: enzyme.get_kcat_values(rxn_ids = [self.rxn_id+"_"+enzyme.id])})
                 continue
 
             self.enzymes.append(enzyme)
@@ -213,23 +229,48 @@ def add_enzymes(self,enzyme_kcat_dict: dict):
             #add constraints to the catalytic event
             self.constraints = {**self.constraints, **enzyme._constraints}
 
+            #add reaction-protein association
+            self.add_enzyme_reaction_association(enzyme)
+
             #connect the enzyme variable to the enzyme in the model and the reaction
             for direction, kcatvalue in kcat.items():
-                coeff = kcatvalue * 3600 * 1e-6
                 #add enzyme to the associated reaction with kinetic constants
                 #and relate enzyme to the catalytic event
                 if direction == 'f':
                     self.constraints[f'EC_{enzyme.id}_{direction}'].set_linear_coefficients({
-                        self.rxn.forward_variable: 1/coeff,
                         enzyme_var.forward_variable: -1
                         })
 
                 elif direction == 'b':
                     self.constraints[f'EC_{enzyme.id}_{direction}'].set_linear_coefficients({
-                        self.rxn.reverse_variable: 1/coeff,
                         enzyme_var.reverse_variable: -1
                         })
 
+    def add_enzyme_reaction_association(self, enzyme):
+        catalytic_reaction = Reaction(self.id + "_" + enzyme.id)
+        self._model.add_reactions([catalytic_reaction])
+        self.catalytic_reactions.append(catalytic_reaction)
+        self._model.constraints[self.id].set_linear_coefficients({
+            catalytic_reaction.forward_variable: -1,
+            catalytic_reaction.reverse_variable: 1,
+        })
+        self.add_catalytic_reaction_to_enzyme_constraint(catalytic_reaction, enzyme)
+
+    def add_catalytic_reaction_to_enzyme_constraint(self, catalytic_reaction:Reaction,
+                                                    enzyme):
+        #remove relation to metabolic reaction
+        for dir in ['f', 'b']:
+            self._model._change_kcat_in_enzyme_constraint(self.rxn, enzyme.id,
+                                                      dir, 0)
+        kcat_dict = self.kcats[enzyme]
+        for direction, kcatvalue in kcat_dict.items():
+            self._model._change_kcat_in_enzyme_constraint(catalytic_reaction, enzyme.id,
+                                                          direction, kcatvalue)
+
+        #change rxn2kcat dict for correct referencing
+        del enzyme.rxn2kcat[self.rxn_id]
+        enzyme.rxn2kcat[catalytic_reaction.id] = kcat_dict
+
 
     def remove_enzymes(self, enzyme_list: list):
         """
@@ -268,18 +309,6 @@ def remove_enzymes(self, enzyme_list: list):
                 self.enzyme_variables.remove(enzyme_var)
             #set coefficient in constraint to 0
             for constraint in set([cons for name, cons in self.constraints.items() if enz.id in name]):
-                # self.constraints[constraint.name] = constraint
-                # coeff = 0
-                # #set coefficients to 0
-                # if constraint.name[-1] == 'f':
-                #     constraint.set_linear_coefficients({
-                #         self.rxn.forward_variable: coeff
-                #         })
-                #
-                # elif constraint.name[-1] == 'b':
-                #     constraint.set_linear_coefficients({
-                #         self.rxn.reverse_variable: coeff
-                #         })
                 if constraint in self._model.constraints.values():
                     self._model.remove_cons_vars([constraint])
                 #remove constraint from list of r=constraints
@@ -330,25 +359,12 @@ def change_kcat_values(self, enzyme_kcat_dict : dict):
 
             enzyme_var = self._model.enzyme_variables.get_by_id(enzyme)
             if enzyme_var not in self.enzyme_variables: self.enzyme_variables.add(enzyme_var)
-            enzyme_obj = self._model.enzymes.get_by_id(enzyme)
-            enzyme_var.change_kcat_values({self.rxn: kcat_dict})
+            ce_reaction = self._model.reactions.get_by_id(self.rxn_id+"_"+enzyme.id)
+            enzyme_var.change_kcat_values({ce_reaction: kcat_dict})
             for direction, kcat in kcats_change.items():
                 #change enzyme variable
-                enzyme_var.kcats[self.rxn_id][direction] = kcat
-                # get constraint
-                constraint_id = f'EC_{enzyme}_{direction}'
-                constraint = enzyme_obj._constraints[constraint_id]
-                # change kcat value in the constraint
-                coeff =  kcat * 3600 * 1e-6
-                if direction == 'f':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        self.rxn.forward_variable: 1/coeff
-                        })
-                elif direction == 'b':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        self.rxn.reverse_variable: 1/coeff
-                        })
-            self._model.solver.update()
+                enzyme_var.kcats[ce_reaction.id][direction] = kcat
+                self._model._change_kcat_in_enzyme_constraint(ce_reaction, enzyme,direction, kcat)
 
     def __copy__(self) -> 'CatalyticEvent':
         """ Copy the CatalyticEvent
diff --git a/src/PAModelpy/Enzyme.py b/src/PAModelpy/Enzyme.py
index 29d67f9..71e40a1 100644
--- a/src/PAModelpy/Enzyme.py
+++ b/src/PAModelpy/Enzyme.py
@@ -4,7 +4,7 @@
 """
 import cobra.core
 import cobra
-from cobra import Reaction, Object
+from cobra import Reaction, Object, Gene, DictList
 from cobra.exceptions import OptimizationError
 from cobra.util.solver import check_solver_status
 from copy import copy, deepcopy
@@ -20,34 +20,23 @@
 
 
 class Enzyme(Object):
-    """Upper level Enzyme object containing information about the enzyme
-    and link to the EnzymeVariables for each reaction the
-    enzyme catalyzes.
-    This class is used to generate enzyme instances from kcat values and contains
-    information about the forward as well as the backward catalysis.
-
-    The enzyme is linked to individual cobra.Reaction variables with CatalyticEvent objects.
-
-    There are two scenarios:
-    - Promiscuous enzymes: a single enzyme can catalyze multiple reactions
-    - Other: a single enzyme catalyzes a single reaction
-
-    Parameters
-    -------
-    id : str
-        Identifier for the enzyme (e.g. Uniprot ID)
-    rxn2kcat: Dict
-        Dictionary with reaction ID, kcat value pairs for the forward (f)
-        and backward (b) reaction
-        (Example: {'PGI': {'f': 30, 'b': 0.1}})
-    upper_bound: float
-        Upper bound for the enzyme variable (default 1000.0)
-    lower_bound: float
-        Lower bound for the enzyme variable (default 0)
-    name: str
-        Name of the enzyme (default None)
-    molmass: float
-        Molar mass of the enzyme (default 3.947778784340140e04)
+    """Upper level Enzyme object containing information about the enzyme and links to the EnzymeVariables for each reaction the enzyme catalyzes.
+
+    Args:
+        - id (str): Identifier for the enzyme (e.g., Uniprot ID).
+        - rxn2kcat (Dict): Dictionary with reaction ID, kcat value pairs for the forward (f) and backward (b) reaction,
+            e.g. `{'PGI': {'f': 30, 'b': 0.1}}`
+        - upper_bound (float): Upper bound for the enzyme variable (default 1000.0).
+        - lower_bound (float): Lower bound for the enzyme variable (default 0).
+        - name (str): Name of the enzyme (default None).
+        - molmass (float): Molar mass of the enzyme (default 3.947778784340140e04).
+
+    Notes:
+        - This class is used to generate enzyme instances from kcat values and contains information about the forward as well as the backward catalysis.
+        - The enzyme is linked to individual cobra.Reaction variables with CatalyticEvent objects.
+        - There are two scenarios:
+            - Promiscuous enzymes: a single enzyme can catalyze multiple reactions.
+            - Other: a single enzyme catalyzes a single reaction.
     """
 
     
@@ -58,6 +47,8 @@ def __init__(
         self,
         id: str,
         rxn2kcat: Dict,
+        rxn2pr: Dict = {},
+        genes: list = [],
         upper_bound: Union[int, float] = 1000.0,
         lower_bound: Union[int, float] = 0,
         name: Optional[str] = None,
@@ -65,6 +56,7 @@ def __init__(
     ):
 
         self.rxn2kcat = rxn2kcat
+        self.rxn2pr = rxn2pr
         self.molmass = molmass
         self.id = id # use e.g. Uniprot ID
         self.name = name
@@ -84,36 +76,35 @@ def __init__(
         self._constraints = {} # dict with constraint_id:optlang.Constraint, key:value pairs.
         self._model = None
         self.enzyme_complex = [] #is the enzyme in a complex?
+        self.genes = genes
+        self.transcripts = DictList()
         self.annotation = {'type':'Constraint'}#you can add an annotation for an enzyme
      
     @property
     def kcat_values(self):
-        """returns a dictionary with kcat values for each associated reaction
+        """Returns a dictionary with kcat values for each associated reaction.
+
+        Returns:
+            dict: A dictionary containing kcat values for associated reactions.
         """
         return self.get_kcat_values()
 
     @property
-    def concentration(self, units:str = 'mmol/gDW', return_units:bool = False) -> float:
-        """returns the enzyme's total concentration
-        Any associated reaction is considered
-
-        Parameters
-        -------
-        units: optional, string
-            units in which the concentration is calculated (default is mmol/gDW), other option is 'g/gDW'
-        return_units: optional, bool
-            determines wheter the units should be returned as well
-
-        Returns
-        -------
-        float
-            Enzyme concentration
+    def concentration(
+        self, units: str = "mmol/gDW", return_units: bool = False
+    ) -> float:
+        """Returns the enzyme's total concentration considering any associated reactions.
+
+        Args:
+            units (str, optional): Units in which the concentration is calculated (default is 'mmol/gDW'), other option is 'g/gDW'.
+            return_units (bool, optional): Determines whether the units should be returned as well.
+
+        Returns:
+            float: Enzyme concentration as a float.
         """
 
         # sum up concentrations (aka fluxes) of all enzyme objects
-        concentration = 0.0
-        for catalytic_event in self.catalytic_events:
-            concentration += catalytic_event.flux()
+        concentration = self.enzyme_variable.concentration
         if units == 'g/gDW':
             #converting mmol to grams of protein:
             # [g] = [mmol]* 1e-3 [mol/mmol] * MW[g/mol]
@@ -122,63 +113,110 @@ def concentration(self, units:str = 'mmol/gDW', return_units:bool = False) -> fl
             return concentration, units
         return concentration
 
+    @concentration.setter
+    def concentration(self, conc):
+        self.enzyme_variable.concentration = conc
+
+    @property
+    def model(self):
+        return self._model
+
+    @model.setter
+    def model(self, model):
+        self._model = model
+
+    def set_forward_concentration(self, conc):
+        self.enzyme_variable.set_forward_concentration(conc)
+
+    def set_reverse_concentration(self, conc):
+        self.enzyme_variable.set_reverse_concentration(conc)
+
     def create_constraint(self, extension: str = None):
-        if extension is None: extension = ''
-        else: extension = '_' + extension
+        if extension is None:
+            extension = ""
+        else:
+            extension = "_" + extension
         # return cobra.core.Metabolite(id = f'EC_{self.id}{extension}', compartment='Enzymes')
-        return optlang.Constraint(Zero, name= f'EC_{self.id}{extension}', lb=0, ub=0)
+        return optlang.Constraint(Zero, name=f"EC_{self.id}{extension}", lb=0, ub=0)
 
     def add_catalytic_event(self, ce: CatalyticEvent, kcats: Dict):
-        """
-        Adding catalytic event associated to a reaction to an enzyme
-        Parameters
-        ----------
-        ce: PAModelpy.Variables.CatalyticEvent
-            The catalytic event object to which the enzyme should be added
-        kcats: dict
-            A list with dicts containing direction, kcat key value pairs
-        Returns
-        -------
+        """Adds a catalytic event associated with a reaction to an enzyme.
+
+        Args:
+            ce (PAModelpy.Variables.CatalyticEvent): A CatalyticEvent object to which the enzyme should be added.
+            kcats (dict): A dictionary containing direction and kcat key-value pairs.
 
+        Returns:
+            NoneType: None
         """
         self.catalytic_events += [ce]
-        self.enzyme_variable.add_catalytic_events([ce],[kcats])
+        self.enzyme_variable.add_catalytic_events([ce], [kcats])
+
+    def add_genes(self, gene_list: list, gene_length:list, relation:str = 'OR') -> None:
+        """
+            Add genes to the enzyme and the model related to the enzyme if applicable
+
+            Args:
+                gene_list (list): A list of gene identifiers to be added.
+                gene_length (list): A list of lengths corresponding to each gene.
+                relation (str, optional): The relationship between genes in gene_list.
+                    Defaults to 'OR'. Possible values: 'OR' or 'AND'.
+
+            Raises:
+                ValueError: If an invalid relation is provided.
+
+            Note:
+                If relation is 'OR', each gene in gene_list will be treated as coding for an individual isozyme
+                If relation is 'AND', all genes in gene_list will be treated as coding for peptides in an enzyme complex
+
+            """
+        # check/correct type of arguments
+        if isinstance(gene_list, str): gene_list = [gene_list]
+        if not isinstance(gene_length, list): gene_length = [gene_length]
+
+        if relation == 'OR':
+            genes_to_add = [[Gene(gene)] for gene in gene_list]
+        elif relation == 'AND':
+            genes_to_add = [[Gene(gene) for gene in gene_list]]
+        else:
+            raise ValueError("Invalid relation. Supported values are 'OR' or 'AND'.")
+        self.genes += genes_to_add
+
+        if self._model is not None:
+            print(genes_to_add)
+            self._model.add_genes(genes = genes_to_add, enzymes = [self], gene_lengths=gene_length)
+
 
     def create_catalytic_event(self, rxn_id: str, kcats: Dict):
-        """creates enzyme variables that link to reactions
-
-        Parameters
-        ----------
-        rxn_id : str
-            ID of the associated reaction in the model
-        kcats : Dict
-            kcat values for the forward and backward reaction
-        Returns
-        -------
-        Variables.CatalyticEvent
-            Enzyme variable object
+        """Creates enzyme variables that link to reactions.
+
+        Args:
+            rxn_id (str): ID of the associated reaction in the model.
+            kcats (Dict): A dictionary containing kcat values for the forward and backward reaction.
+
+        Returns:
+            Variables.CatalyticEvent: Enzyme variable object of type Variables.CatalyticEvent.
         """
         
         # create unique enzyme object name and id
         catalytic_event_id = self.catalytic_event_id.format(rxn_id)
         if self.name is not None:
-            enzyme_object_name = rxn_id + '_' + self.name 
+            enzyme_object_name = rxn_id + "_" + self.name
         else:
-            enzyme_object_name = rxn_id + '_'  + self.id
-        
+            enzyme_object_name = rxn_id + "_" + self.id
+
         # create enzymatic reaction object inherited from cobra.Reaction
-        catalytic_event= CatalyticEvent(
+        catalytic_event = CatalyticEvent(
             id=catalytic_event_id,
             rxn_id=rxn_id,
             kcats2enzymes={self: kcats},
-            name=enzyme_object_name
+            name=enzyme_object_name,
         )
 
         self.add_catalytic_event(catalytic_event, kcats)
 
     def create_enzyme_variable(self):
-        """creates enzyme variables that link  enzyme to reactions
-        """
+        """Creates enzyme variables that link  enzyme to reactions."""
         # create enzymatic reaction object inherited from cobra.Reaction
         enzyme_variable = EnzymeVariable(
             id=self.id,
@@ -189,43 +227,35 @@ def create_enzyme_variable(self):
         self.enzyme_variable = enzyme_variable
 
     def change_kcat_values(self, rxn2kcat: Dict):
-        """changes the kcat values for the enzyme
-        and updates the enzyme variable (enzymatic reaction) accordingly
-
-        Parameters
-        ----------
-        rxn2kcat : Dict
-            Dictionary with reaction ID, kcat value pairs for the forward (f)
-            and backward (b) reaction
-            (Example: {'PGI': {'f': 30, 'b': 0.1}})
+        """Changes the kcat values for the enzyme and updates the enzyme variable (enzymatic reaction) accordingly.
+
+        Args:
+            rxn2kcat (Dict): A dictionary with reaction ID, kcat value pairs for the forward (f) and backward (b) reaction,
+                e.g. `{'PGI': {'f': 30, 'b': 0.1}}`
         """
 
         # update the enzyme variables
         for rxn_id, kcats in rxn2kcat.items():
             catalytic_event_id = self.catalytic_event_id.format(rxn_id)
-            #change rxn2kcat dictionary
+            # change rxn2kcat dictionary
             self.rxn2kcat[rxn_id] = kcats
             # is there already a link between enzyme and reaction?
             if catalytic_event_id not in self.catalytic_events:
-                warn(f'Reaction {rxn_id} is not associated with enzyme {self.id}. Skip')
+                warn(f"Reaction {rxn_id} is not associated with enzyme {self.id}. Skip")
             else:
                 # change kcat values of existing enzyme variable
-                self.catalytic_events.get_by_id(catalytic_event_id).change_kcat_values({self.id: kcats})
-
+                self.catalytic_events.get_by_id(catalytic_event_id).change_kcat_values(
+                    {self.id: kcats}
+                )
 
     def get_kcat_values(self, rxn_ids: Union[str, list] = None) -> Dict:
-        """returns the kcat values for a specific enzyme and all
-        enzyme-associated reactions
-
-        Parameters
-        ----------
-        rxn_ids : str or list
-            ID of the reactions for which the kcat values should be returned
-
-        Returns
-        -------
-        Dict
-            kcat values for the forward and backward reaction
+        """Returns the kcat values for a specific enzyme and all enzyme-associated reactions.
+
+        Args:
+            rxn_ids (Union[str, list], optional): ID of the reactions for which the kcat values should be returned. It can be a single reaction ID (str) or a list of reaction IDs.
+
+        Returns:
+            Dict: A dictionary containing kcat values for the forward (f) and backward (b) reactions.
         """
         
         if isinstance(rxn_ids, str):
@@ -240,50 +270,51 @@ def get_kcat_values(self, rxn_ids: Union[str, list] = None) -> Dict:
             for rxn_id in rxn_ids:
                 catalytic_event_id = self.catalytic_event_id.format(rxn_id)
                 if catalytic_event_id not in self.catalytic_events:
-                    warn('No reaction {0} found'.format(rxn_id))
+                    warn("No reaction {0} found".format(rxn_id))
                 else:
                     # get kcat values
                     rxn2kcat[rxn_id] = self.rxn2kcat[rxn_id]
 
-        if len(rxn_ids)==1:
+        if len(rxn_ids) == 1:
             return rxn2kcat[rxn_id]
                     
         return rxn2kcat
 
     def remove_catalytic_event(self, catalytic_event: Union[CatalyticEvent, str]):
-        """
-        Function to remove a catalytic event from an enzyme
-        Parameters
-        ----------
-        catalytic_event: CatalyticEvent or str
-            catalytic event or identifier to remove
+        """Function to remove a catalytic event from an enzyme.
+
+        Args:
+            catalytic_event (Union[CatalyticEvent, str]): CatalyticEvent or str, catalytic event or identifier to remove.
         """
         if isinstance(catalytic_event, str):
             try:
                 catalytic_event = self.catalytic_events.get_by_id(catalytic_event)
             except:
-                print(f'Catalytic event {catalytic_event} is not related to this enzyme and can thus not be removed!')
+                print(
+                    f"Catalytic event {catalytic_event} is not related to this enzyme and can thus not be removed!"
+                )
 
-        #remove the event from the DictList
+        # remove the event from the DictList
         self.catalytic_events.remove(catalytic_event)
 
-    def __copy__(self) -> 'Enzyme':
-        """ Copy the enzyme variable
-        :return: PAModelpy.Enzyme.Enzyme:
-        A new enzyme that is a copy of the original enzyme
+    def __copy__(self) -> "Enzyme":
+        """Copy the enzyme variable.
+
+        Returns:
+            PAModelpy.Enzyme.Enzyme: A new enzyme that is a copy of the original enzyme.
         """
 
         cop = copy(super(Enzyme, self))
         return cop
 
-    def __deepcopy__(self, memo: dict) -> 'Enzyme':
-        """ Copy the enzyme variable with memo
+    def __deepcopy__(self, memo: dict) -> "Enzyme":
+        """Copy the enzyme variable with memo.
 
-        :param: memo:dict:
-        Automatically passed parameter
+        Args:
+            memo (dict): Automatically passed parameter.
 
-        :return: PAModelpy.Enzyme.Enzyme:
-        A new enzyme that is a copy of the original enzyme with memo
+        Returns:
+            PAModelpy.Enzyme.Enzyme: A new enzyme that is a copy of the original enzyme with memo.
         """
 
         cop = deepcopy(super(Enzyme, self), memo)
@@ -291,34 +322,21 @@ def __deepcopy__(self, memo: dict) -> 'Enzyme':
 
 
 class EnzymeComplex(Enzyme):
-    """Upper level EnzymeComplex object containing information about the enzymes in a complex
-       and link to the enzyme variables (CatalyticEvents) for each reaction the
-       enzyme complex catalyzes.
-       This class is used to generate enzyme instances from kcat values and contains
-       information about the forward as well as the backward catalysis.
-
-       There are two scenarios:
-       - Promiscuous enzymes: a single enzyme complex can catalyze multiple reactions
-       - Other: a single enzyme complex catalyzes a single reaction
-
-       Parameters
-       -------
-       id : str
-           Identifier for the enzyme complex (e.g. Uniprot ID)
-       enzymes: DictList of cobra.core.Enzyme
-            Enzyme objects associated with the enzyme complex
-       rxn2kcat: Dict
-           Dictionary with reaction ID, kcat value pairs for the forward (f)
-           and backward (b) reaction
-           (Example: {'PGI': {'f': 30, 'b': 0.1}})
-       upper_bound: float
-           Upper bound for the enzyme variable (default 1000.0)
-       name: str
-           Name of the enzyme (default None)
-       molmass: float
-           Molar mass of the enzyme (default 3.947778784340140e04)
-
-       """
+    """Upper-level EnzymeComplex object containing information about the enzymes in a complex
+    and a link to the enzyme variables (CatalyticEvents) for each reaction the enzyme complex catalyzes.
+
+    This class is used to generate enzyme instances from kcat values and contains
+    information about the forward as well as the backward catalysis.
+
+    Args:
+        id (str): Identifier for the enzyme complex (e.g., Uniprot ID).
+        enzymes (DictList of cobra.core.Enzyme): Enzyme objects associated with the enzyme complex.
+        rxn2kcat (Dict): Dictionary with reaction ID, kcat value pairs for the forward (f) and backward (b) reaction,
+            e.g. `{'PGI': {'f': 30, 'b': 0.1}}`
+        upper_bound (float, optional): Upper bound for the enzyme variable (default 1000.0).
+        name (str, optional): Name of the enzyme (default None).
+        molmass (float, optional): Molar mass of the enzyme (default 3.947778784340140e04).
+    """
 
     # constant parameters
     DEFAULT_ENZYME_MOL_MASS = 3.947778784340140e04  # mean enzymes mass E.coli [g/mol]
@@ -326,59 +344,59 @@ class EnzymeComplex(Enzyme):
     def __init__(
             self,
             id: str,
-            enzymes: DictList,
             rxn2kcat: Dict,
+            enzymes: DictList = DictList(),
+            genes: list = [],
             upper_bound: Union[int, float] = 1000.0,
             name: Optional[str] = None,
-            molmass: Union[int, float] = DEFAULT_ENZYME_MOL_MASS,
-    ):
+            molmass: Union[int, float] = DEFAULT_ENZYME_MOL_MASS, ):
         super().__init__(
             id = id,
             rxn2kcat=rxn2kcat,
-            upper_bound = upper_bound,
+            genes=genes,
+            upper_bound=upper_bound,
             name=name,
             molmass=molmass
         )
 
-        self.enzymes = None
+        self.reactions = list(rxn2kcat.keys())
+        self.enzymes = DictList()
         self.add_enzymes(enzymes)
 
     def add_enzymes(self, enzymes: DictList):
         for enzyme in enzymes:
-            self.enzymes.append(enzyme)
-            enzyme.enzyme_complex.append(self.id)
-            self.molmass += enzyme.molmass
-
+            if enzyme not in self.enzymes:
+                self.enzymes.append(enzyme)
+                enzyme.enzyme_complex.append(self.id)
+                # self.molmass += enzyme.molmass
+                #remove associated constraints from the model
+                if self._model is not None:
+                    if enzyme in self._model.enzymes:
+                        for rxn in self.reactions:
+                            self._model.remove_enzyme_reaction_association(enzyme, rxn)
+                    else:
+                        self._model.enzymes.append(enzyme)
 
 class EnzymeVariable(Reaction):
+    """EnzymeVariable is a class for holding information regarding the variable representing an enzyme in the model.
+    For each reaction, the enzyme variables are summarized in a CatalyticEvent.
+
+    There are three different scenarios:
+    - Enzyme complex: multiple enzymes together are associated with an EnzymeComplex object.
+    - Isozymes: multiple enzymes independently associated with a single catalytic event.
+    - Other: a single enzyme is associated with a single catalytic event.
+
+    Args:
+        kcats2rxns (Dict): A dictionary with reaction_id, kcat key, value pairs to connect the enzyme with the associated reaction.
+            The kcat is another dictionary with `f` and `b` for the forward and backward reactions, respectively.
+        id (str, optional): The identifier to associate with this enzyme (default None).
+        name (str, optional): A human-readable name for the enzyme (default "").
+        subsystem (str, optional): Subsystem where the enzyme is meant to function (default "").
+        lower_bound (float): The lower flux bound (default 0.0).
+        upper_bound (float, optional): The upper flux bound (default None).
+        **kwargs: Additional keyword arguments are passed on to the parent class.
     """
-           EnzymeVariable is a class for holding information regarding the
-           variable representing an enzyme in the model. For each reaction, the enzyme variables are
-           summarized in a CatalyticEvent
-           There are three different scenarios:
-           - Enzyme complex: multiple enzymes together are associated with an EnzymeComplex object
-           - isozymes: multiple enzymes independently associated with a single catalytic event
-           - Other: a single enzyme is associated with a single catalytic event
-
-           Parameters
-           ----------
-           kcats2rxns: Dict
-               A Dict with reaction_id, kcat key, value pairs to connect the
-               enzyme with the associated reaction the kcat is another dict with 'f' and 'b'
-               for the forward and backward reactions respectively.
-           id : str, optional
-               The identifier to associate with this enzyme (default None)
-           name : str, optional
-               A human-readable name for the reaction (default "").
-           subsystem : str, optional
-               Subsystem where the reaction is meant to occur (default "").
-           lower_bound : float
-               The lower flux bound (default 0.0).
-           upper_bound : float, optional
-               The upper flux bound (default None).
-           **kwargs:
-               Further keyword arguments are passed on to the parent class.
-           """
+
     DEFAULT_ENZYME_MOL_MASS = 3.947778784340140e04  # mean enzymes mass E.coli [g/mol]
 
     def __init__(
@@ -401,59 +419,56 @@ def __init__(
         self.kcats = kcats2rxns
         self.molmass = molmass
         self.rxn_ids = [rxn for rxn in kcats2rxns.keys()]
-        self.enzyme = None #store the enzyme associated to the enzyme variable
-        self.catalytic_events = DictList() #store the interfaces to reactions related to this enzyme
+        self.enzyme = None  # store the enzyme associated to the enzyme variable
+        self.catalytic_events = (
+            DictList()
+        )  # store the interfaces to reactions related to this enzyme
         self.reactions = DictList()
         self.constraints = {}  # store IDs of constraint the enzyme is associated with
         self._model = None
-        self.annotation = {'type': 'Variable'}
+        self.annotation = {"type": "Variable"}
         self.variables = dict()
 
     @property
     def kcat_values(self):
-        """returns a dictionary with kcat values and reactions
+        """Returns a dictionary with kcat values and reactions.
+
+        Returns:
+            dict: A dictionary containing kcat values and their associated reactions.
         """
         return self.kcats
 
     @property
     def flux(self) -> float:
-        """
-        Get the flux value in the most recent solution.
+        """Get the flux value in the most recent solution.
+
         Flux is the primal value of the corresponding variable in the model.
-        Returns
-        -------
-        flux: float
-            Flux is the primal value of the corresponding variable in the model.
-        Warnings
-        --------
-        * Accessing reaction fluxes through a `Solution` object is the safer,
-          preferred, and only guaranteed to be correct way. You can see how to
-          do so easily in the examples.
-        * Reaction flux is retrieved from the currently defined
-          `self._model.solver`. The solver status is checked but there are no
-          guarantees that the current solver state is the one you are looking
-          for.
-        * If you modify the underlying model after an optimization, you will
-          retrieve the old optimization values.
-        Raises
-        ------
-        RuntimeError
-            If the underlying model was never optimized beforehand or the
-            reaction is not part of a model.
-        OptimizationError
-            If the solver status is anything other than 'optimal'.
-        AssertionError
-            If the flux value is not within the bounds.
-        Examples
-        --------
-        >>> from cobra.io import load_model
-        >>> model = load_model("textbook")
-        >>> solution = model.optimize()
-        >>> model.variables.PFK.flux
-        7.477381962160283
-        >>> solution.fluxes.PFK
-        7.4773819621602833
+
+        Returns:
+            float: Flux, which is the primal value of the corresponding variable in the model.
+
+        Raises:
+            RuntimeError: If the underlying model was never optimized beforehand or the reaction is not part of a model.
+            OptimizationError: If the solver status is anything other than 'optimal'.
+            AssertionError: If the flux value is not within the bounds.
+
+        Warnings:
+            - Accessing reaction fluxes through a `Solution` object is the safer, preferred, and only guaranteed to be correct way.
+            - Reaction flux is retrieved from the currently defined `self._model.solver`. The solver status is checked but there are no guarantees that the current solver state is the one you are looking for.
+            - If you modify the underlying model after an optimization, you will retrieve the old optimization values.
+
+        Examples:
+            ```
+            >>> from cobra.io import load_model
+            >>> model = load_model("textbook")
+            >>> solution = model.optimize()
+            >>> model.variables.PFK.flux
+            7.477381962160283
+            >>> solution.fluxes.PFK
+            7.4773819621602833
+            ```
         """
+
         try:
             check_solver_status(self._model.solver.status)
             return self.forward_variable.primal + self.reverse_variable.primal
@@ -470,16 +485,59 @@ def flux(self) -> float:
 
     @property
     def concentration(self) -> float:
+        """Get the enzyme concentration value of the most recent solution.
+
+        The enzyme concentration equals the flux value.
+
+        Returns:
+            float: Enzyme concentration in mmol/gDW."""
+        return self.flux*1e-6
+
+    @concentration.setter
+    def concentration(self, conc:Union[float,int]) -> None:
         """
-        Get the enzyme concentration value of the most recent solution.
-        The enzyme concentration equals the flux value
+        Sets the concentration of the enzyme by creating or updating a constraint
+        that enforces the concentration to be equal to the sum of the forward and reverse
+        variable primals.
+
+        Args:
+        conc : float, int
+            The concentration value to be set for the enzyme. This value will be used
+            as both the lower and upper bound for the constraint, effectively fixing the
+            concentration to this value.
+
+        Notes
+        -----
+        - If a concentration constraint for the enzyme does not already exist in the model,
+          this function creates a new constraint named '<enzyme_id>_conc'.
+        - The concentration constraint is defined as:
+          concentration = forward_variable.primal + reverse_variable.primal
+        - If the constraint already exists, the linear coefficients for the forward and reverse
+          variables are updated to ensure the constraint remains valid.
 
-        Returns
-        -------
-        float
-            enzyme concentration [mmol/gDW]
+        Raises
+        ------
+        ValueError
+            If `conc` is not a valid numerical value.
         """
-        return self.flux
+
+        # Ensure `conc` is a valid numerical value
+        if not isinstance(conc, (int, float)):
+            raise ValueError("The concentration must be a numerical value.")
+
+        # Make the constraint: concentration = forward_variable + reverse_variable
+        if self.id + '_conc' not in self._model.constraints.keys():
+            concentration_constraint = self._model.problem.Constraint(
+                Zero, name=self.id + '_conc', lb=conc, ub=conc)
+            self._model.add_cons_vars(concentration_constraint)
+
+        #units of enzyme should be *1e6 because of feasibility tolerance solver
+        self._model.constraints[self.id + '_conc'].set_linear_coefficients({
+            self.forward_variable: 1e-6,
+            self.reverse_variable: 1e-6
+        })
+        self.enzyme._constraints[self.id+'_conc'] = self._model.constraints[self.id + '_conc']
+
 
     @property
     def model(self):
@@ -502,9 +560,16 @@ def model(self, model):
         self._model.enzyme_variables.append(self)
 
         # create new forward and reverse variables
-        forward_variable = self._model.problem.Variable(name=self.id, ub=self.upper_bound, lb=0)
-        reverse_variable = self._model.problem.Variable(name=self.reverse_id, ub=-self.lower_bound, lb=0)
-        self.variables = {'forward_variable': forward_variable, 'reverse_variable': reverse_variable}
+        forward_variable = self._model.problem.Variable(
+            name=self.id, ub=self.upper_bound, lb=0
+        )
+        reverse_variable = self._model.problem.Variable(
+            name=self.reverse_id, ub=-self.lower_bound, lb=0
+        )
+        self.variables = {
+            "forward_variable": forward_variable,
+            "reverse_variable": reverse_variable,
+        }
         self._model.add_cons_vars([forward_variable, reverse_variable])
 
 
@@ -516,51 +581,118 @@ def model(self, model):
                 self._model.add_reactions([rxn])
             if not rxn_id in self.reactions:
                 self.reactions.append(self._model.reactions.get_by_id(rxn_id))
-            #create a catalytic event if it doesn't exist already
-            if not f'CE_{rxn_id}' in self._model.catalytic_events:
+            # create a catalytic event if it doesn't exist already
+            if not f"CE_{rxn_id}" in self._model.catalytic_events:
                 kcats2enzymes = {self.enzyme: self.kcats[rxn_id]}
-                ce = CatalyticEvent(id = f'CE_{rxn_id}',
-                                    kcats2enzymes=kcats2enzymes,
-                                    rxn_id=rxn_id)
+                ce = CatalyticEvent(
+                    id=f"CE_{rxn_id}", kcats2enzymes=kcats2enzymes, rxn_id=rxn_id
+                )
                 ce.model = model
                 self._model.catalytic_events.append(ce)
             # if catalytic event exist, add the enzyme to it
             else:
-                self._model.catalytic_events.get_by_id(f'CE_{rxn_id}').add_enzymes({self.enzyme: self.kcats[rxn_id]})
+                self._model.catalytic_events.get_by_id(f"CE_{rxn_id}").add_enzymes(
+                    {self.enzyme: self.kcats[rxn_id]}
+                )
             # if catalytic event is not related to the enzyme variable yet, add it.
-            if f'CE_{rxn_id}' not in self.catalytic_events:
-                self.catalytic_events.append(self._model.catalytic_events.get_by_id(f'CE_{rxn_id}'))
+            if f"CE_{rxn_id}" not in self.catalytic_events:
+                self.catalytic_events.append(
+                    self._model.catalytic_events.get_by_id(f"CE_{rxn_id}")
+                )
 
     @property
     def forward_variable(self):
         if self._model is not None:
             return self._model.variables[self.id]
         else:
-            return self.variables['forward_variable']
+            return None
 
     @property
     def reverse_variable(self):
         if self._model is not None:
             return self._model.variables[self.reverse_id]
         else:
-            return self.variables['reverse_variable']
+            return None
+
+    def set_forward_concentration(self, conc: Union[float, int]) -> None:
+        """
+        Sets the concentration of the enzyme by creating or updating a constraint
+        that enforces the concentration to be equal to the sum of only the forward
+        variable primals. This forces a reaction to be active in the forward direction.
+        It used the concentration setter functionality and subsequently sets the
+        coefficient for the reverse variable in the constraint to 0.
+
+        Args:
+        conc : float, int
+            The concentration value to be set for the enzyme. This value will be used
+            as both the lower and upper bound for the constraint, effectively fixing the
+            concentration to this value.
+
+        Notes
+        -----
+        - If a concentration constraint for the enzyme does not already exist in the model,
+          this function creates a new constraint named '<enzyme_id>_conc'.
+        - The concentration constraint is defined as:
+          concentration = forward_variable.primal
+
+        Raises
+        ------
+        ValueError
+            If `conc` is not a valid numerical value.
+        """
+        self.concentration = conc
+        self._model.constraints[self.id + '_conc'].set_linear_coefficients({
+            self.forward_variable: 1,
+            self.reverse_variable: 0
+        })
+
+    def set_reverse_concentration(self, conc: Union[float, int]) -> None:
+        """
+        Sets the concentration of the enzyme by creating or updating a constraint
+        that enforces the concentration to be equal to the sum of only the reverse
+        variable primals. This forces a reaction to be active in the reverse direction.
+        It used the concentration setter functionality and subsequently sets the
+        coefficient for the forward variable in the constraint to 0.
+
+        Args:
+        conc : float, int
+            The concentration value to be set for the enzyme. This value will be used
+            as both the lower and upper bound for the constraint, effectively fixing the
+            concentration to this value.
+
+        Notes
+        -----
+        - If a concentration constraint for the enzyme does not already exist in the model,
+          this function creates a new constraint named '<enzyme_id>_conc'.
+        - The concentration constraint is defined as:
+          concentration = reverse_variable.primal
+
+        Raises
+        ------
+        ValueError
+            If `conc` is not a valid numerical value.
+        """
+        self.concentration = conc
+        self._model.constraints[self.id + '_conc'].set_linear_coefficients({
+            self.forward_variable: 0,
+            self.reverse_variable: 1
+        })
 
     def add_catalytic_events(self, catalytic_events:list, kcats:list):
         """
         Adding a catalytic event to an enzyme variable
 
-        Parameters
-        ----------
-        catalytic_events: list
-            Catalytic events to add
-        kcats:list
-            A list with dicts containing direction, kcat key value pairs
+        Args:
+            catalytic_events (list): A list of catalytic events to add.
+            kcats (list): A list with dictionaries containing direction and kcat key-value pairs.
         """
 
         for i, ce in enumerate(catalytic_events):
             if ce in self.catalytic_events:
-                warn(f'Catalytic event {ce.id} is already associated with enzyme variable {self.id}. '
-                     f'Continue with other catalytic events')
+                warn(
+                    f"Catalytic event {ce.id} is already associated with enzyme variable {self.id}. "
+                    f"Continue with other catalytic events"
+                )
             else:
                 if not ce.enzyme_variables.has_id(self.id):
                     ce.enzyme_variables.append(self)
@@ -569,25 +701,21 @@ def add_catalytic_events(self, catalytic_events:list, kcats:list):
                 self.add_reactions({ce.rxn_id: kcats[i]})
 
     def add_reactions(self, reaction_kcat_dict: dict):
-        """
-        Add reactions to the enzyme variable and create bindings to the related model.
+        """Add reactions to the enzyme variable and create bindings to the related model.
         If there are multiple reactions related to a single enzyme, this is an isozyme.
 
-        Parameters
-        ----------
-        reaction_kcat_dict: nested dict
-            A Dict with the reaction_id, kcat key, value pairs to connect the
-            enzyme with the associated reaction the kcat is another dict with 'f' and 'b'
-            for the forward and backward reactions respectively.
-
-
+        Args:
+            reaction_kcat_dict (dict): A nested dictionary with the reaction_id, kcat key, value pairs to connect the
+                enzyme with the associated reaction. The kcat is another dictionary with `f` and `b` for the forward and
+                backward reactions, respectively.
         """
         for reaction, kcat in reaction_kcat_dict.items():
             # check if the enzyme is already associated to the catalytic event
             try:
                 self.reactions.get_by_id(reaction.id)
                 warn(
-                    f'Reaction {reaction.id} is already associated with the enzyme {self.id}. The enzyme variable will be updated')
+                    f"Reaction {reaction.id} is already associated with the enzyme {self.id}. The enzyme variable will be updated"
+                )
                 self.change_kcat_values(kcat)
                 return
             except:
@@ -610,54 +738,52 @@ def add_reactions(self, reaction_kcat_dict: dict):
             self.reactions.append(rxn)
 
             for direction in kcat.keys():
-                if direction != 'f' and direction != 'b':
-                    warn(f'Invalid direction {direction} for kcat value for enzyme variable {self.id}! Skip!')
+                if direction != "f" and direction != "b":
+                    warn(
+                        f"Invalid direction {direction} for kcat value for enzyme variable {self.id}! Skip!"
+                    )
                     continue
 
             # add enzyme to catalytic event and the related variable
             for direction, kcatvalue in kcat.items():
-                coeff = kcatvalue * 3600 * 1e-6
+                self._model._change_kcat_in_enzyme_constraint(rxn, self.id,
+                                                              direction, kcatvalue)
                 # add enzyme to the associated reaction with kinetic constants
                 # and relate enzyme to the catalytic event
                 if direction == 'f':
                     self.constraints[f'EC_{self.id}_{direction}'].set_linear_coefficients({
-                        rxn.forward_variable: 1 / coeff,
                         self.forward_variable: -1
                     })
 
                 elif direction == 'b':
                     self.constraints[f'EC_{self.id}_{direction}'].set_linear_coefficients({
-                        rxn.reverse_variable: 1 / coeff,
                         self.reverse_variable: -1
                     })
 
     def remove_catalytic_event(self, catalytic_event: Union[CatalyticEvent, str]):
-        """
-        Function to remove a catalytic event from an enzyme
-        Parameters
-        ----------
-        catalytic_event: CatalyticEvent or str
-            catalytic event or identifier to remove
+        """Remove a catalytic event from an enzyme.
+
+        Args:
+            catalytic_event (Union[CatalyticEvent, str]): CatalyticEvent or str, catalytic event or identifier to remove.
         """
         if isinstance(catalytic_event, str):
             try:
                 catalytic_event = self.catalytic_events.get_by_id(catalytic_event)
             except:
-                print(f'Catalytic event {catalytic_event} is not related to this enzyme and can thus not be removed!')
+                print(
+                    f"Catalytic event {catalytic_event} is not related to this enzyme and can thus not be removed!"
+                )
 
-        #remove the event from the DictList
+        # remove the event from the DictList
         self.catalytic_events.remove(catalytic_event.id)
 
     def remove_reactions(self, reaction_list: list):
-        """
-        Remove reaction from the enzyme variable and remove the reaction from the
-        constraint expressions related to the enzyme
-
-        Parameters
-        ----------
-        reaction_list: list
-            A list with Cbra.Reaction objects which should be removed. If a list of identifiers (str)
-            is provided, the corresponding enzyme will be obtained from the EnzymeVariables.reaction attribute
+        """Remove reactions from the enzyme variable and remove the reactions from the
+        constraint expressions related to the enzyme.
+
+        Args:
+            reaction_list (list): A list with Cobra.Reaction objects which should be removed. If a list of identifiers (str)
+                is provided, the corresponding enzyme will be obtained from the EnzymeVariables.reaction attribute.
         """
         # check the input
         if not hasattr(reaction_list, "__iter__"):
@@ -672,7 +798,8 @@ def remove_reactions(self, reaction_list: list):
                     reaction_list[i] = self.reactions.get_by_id(rxn)
                 except:
                     print(
-                        f'Reaction {rxn} is not associated with the enzyme variable {self.id}. This reaction cannot be removed. \n')
+                        f"Reaction {rxn} is not associated with the enzyme variable {self.id}. This reaction cannot be removed. \n"
+                    )
                     pass
 
         for rxn in reaction_list:
@@ -685,29 +812,25 @@ def remove_reactions(self, reaction_list: list):
                 self.constraints[constraint.name] = constraint
                 coeff = 0
                 # set coefficients to 0
-                if constraint.name[-1] == 'f':
-                    constraint.set_linear_coefficients({
-                        rxn.forward_variable: coeff,
-                        self.forward_variable: 0
-                    })
-
-                elif constraint.name[-1] == 'b':
-                    constraint.set_linear_coefficients({
-                        rxn.reverse_variable: coeff,
-                        self.reverse_variable: 0
-                    })
+                if constraint.name[-1] == "f":
+                    constraint.set_linear_coefficients(
+                        {rxn.forward_variable: coeff, self.forward_variable: 0}
+                    )
+
+                elif constraint.name[-1] == "b":
+                    constraint.set_linear_coefficients(
+                        {rxn.reverse_variable: coeff, self.reverse_variable: 0}
+                    )
                 # remove constraint from list of r=constraints
                 del self.constraints[constraint.name]
 
     def change_kcat_values(self, reaction_kcat_dict: dict):
-        """changes kcat values for the enzyme variable
-        Parameters
-        ----------
-        reaction_kcat_dict: nested Dict
-            A Dict with Cobra.Reaction, kcat key, value pairs to connect the
-            enzyme with the associated reaction the kcat is another dict with 'f' and 'b'
-            for the forward and backward reactions respectively.
+        """Changes kcat values for the enzyme variable.
 
+        Args:
+            reaction_kcat_dict (dict): A nested dictionary with Cobra.Reaction, kcat key, value pairs to connect the
+                enzyme with the associated reaction. The kcat is another dictionary with `f` and `b` for the forward and
+                backward reactions, respectively.
         """
         # apply changes to internal dicts (one by one to avoid deleting kcat values)
         kcats_change = {}
@@ -715,49 +838,41 @@ def change_kcat_values(self, reaction_kcat_dict: dict):
             # save change in dict
             self.kcats[rxn.id] = kcat_dict
             for direction, kcat in kcat_dict.items():
-                if direction != 'f' and direction != 'b':
-                    warn(f'Invalid direction {direction} for kcat value for enzyme variable {self.id}! Skip!')
+                if direction != "f" and direction != "b":
+                    warn(
+                        f"Invalid direction {direction} for kcat value for enzyme variable {self.id}! Skip!"
+                    )
                     continue
 
                 kcats_change[direction] = kcat
 
             # is enzyme variable already integrated into a model
             if self._model is None:
-                warn(f'Catalytic event {self.id} is not integrated into a model!')
+                warn(f"Catalytic event {self.id} is not integrated into a model!")
 
             for direction, kcat in kcats_change.items():
-                # get constraint
-                constraint_id = f'EC_{self.id}_{direction}'
-                constraint = self.enzyme._constraints[constraint_id]
-                # change kcat value in the constraint
-                coeff = kcat * 3600 * 1e-6
-                if direction == 'f':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        rxn.forward_variable: 1 / coeff
-                    })
-                elif direction == 'b':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        rxn.reverse_variable: 1 / coeff
-                    })
+                self._model._change_kcat_in_enzyme_constraint(rxn, self.id,
+                                                              direction, kcat)
             self._model.solver.update()
 
-    def __copy__(self) -> 'PAModelpy.Enzyme.EnzymeVariable':
-        """ Copy the enzyme variable
-        :return: PAModelpy.Enzyme.EnzymeVariable:
-        A new enzyme variable that is a copy of the original enzyme variable
+    def __copy__(self) -> "PAModelpy.Enzyme.EnzymeVariable":
+        """Copy the enzyme variable.
+
+        Returns:
+            PAModelpy.Enzyme.EnzymeVariable: A new enzyme variable that is a copy of the original enzyme variable.
         """
 
         cop = copy(super(EnzymeVariable, self))
         return cop
 
-    def __deepcopy__(self, memo: dict) -> 'PAModelpy.Enzyme.EnzymeVariable':
-        """ Copy the enzyme variable with memo
+    def __deepcopy__(self, memo: dict) -> "PAModelpy.Enzyme.EnzymeVariable":
+        """Copy the enzyme variable with memo.
 
-        :param: memo:dict:
-        Automatically passed parameter
+        Args:
+            memo (dict): Automatically passed parameter.
 
-        :return: PAModelpy.Enzyme.EnzymeVariable:
-        A new enzyme variable that is a copy of the original enzyme variable with memo
+        Returns:
+            PAModelpy.Enzyme.EnzymeVariable: A new enzyme variable that is a copy of the original enzyme variable with memo.
         """
 
         cop = deepcopy(super(EnzymeVariable, self), memo)
diff --git a/src/PAModelpy/EnzymeSectors.py b/src/PAModelpy/EnzymeSectors.py
index d9139b3..2497d49 100644
--- a/src/PAModelpy/EnzymeSectors.py
+++ b/src/PAModelpy/EnzymeSectors.py
@@ -1,8 +1,9 @@
 from warnings import warn
 from copy import copy, deepcopy
-from cobra import Object
+from cobra import Object, Gene, Model
+from typing import Union
 
-from .Enzyme import Enzyme
+from .Enzyme import Enzyme, EnzymeComplex
 from .configuration import Config
 
 
@@ -58,7 +59,7 @@ def __init__(self, id_list, mol_mass, configuration = Config):
 class ActiveEnzymeSector(Sector):
     DEFAULT_MOL_MASS = 3.947778784340140e04  # mean enzymes mass E.coli [g/mol]
     #class with all the information on enzymes related to metabolic reactions
-    def __init__(self, rxn2protein: dict, configuration:Config =Config):
+    def __init__(self, rxn2protein: dict, protein2gene:dict = {},configuration:Config =Config):
         """_summary_
 
         Parameters
@@ -72,10 +73,19 @@ def __init__(self, rxn2protein: dict, configuration:Config =Config):
               example: {R1:
                             {E1:
                                 {'f': forward kcat, 'b': backward kcat, 'molmass': molar mass, 'genes': [G1, G2],
+                                'protein_reaction_association': [[E1, E2]],
                                  'complex_with': 'E2'},
                             E2:
                                 {'f': forward kcat, 'b': backward kcat, 'molmass': molar mass, 'genes': [G3, G4],
-                                 'complex_with': 'E1'}}
+                                 'protein_reaction_association': [[E1, E2]], 'complex_with': 'E1'}}
+        protein2gene: dict
+            enzyme_id, gene_list key, value pairs for each enzyme.The gene_list value is a list of lists which indicates
+            'and' or 'or' relationships between the genes which code for the enzyme(complex).
+
+            example: {E1: [[gene1], [gene2, gene3]], E2: [[gene4]]}
+            where the gene-protein-reaction associations are the following:
+            E1: gene1 or (gene2 and gene3)
+            E2: gene4
 
         configuration: Config object, optional
                     Information about general configuration of the model including identifier conventions.
@@ -84,6 +94,7 @@ def __init__(self, rxn2protein: dict, configuration:Config =Config):
         self.TOTAL_PROTEIN_CONSTRAINT_ID = configuration.TOTAL_PROTEIN_CONSTRAINT_ID
 
         self.rxn2protein = rxn2protein
+        self.protein2gene = protein2gene
         # ID of the sector
         self.id = 'ActiveEnzymeSector'
         self.model = None
@@ -129,6 +140,7 @@ def add(self, model):
 
             for enzyme_id, enzyme_dict in enzymes.items():
                 kcat = dict((k, v) for k, v in enzyme_dict.items() if k == 'f' or k == 'b')
+                protein_reaction = enzyme_dict['protein_reaction_association']
 
                 # get molar mass of enzyme or replace with default value
                 if 'molmass' in enzyme_dict.keys():
@@ -138,36 +150,57 @@ def add(self, model):
 
                 # check if there already exists an Enzyme object for the EC numbers associated to this reaction,
                 # otherwise create one and store all information
-                if enzyme_id in model.enzymes:
+                if enzyme_id in model.enzyme_variables and not self._enzyme_is_enzyme_complex(protein_reaction, enzyme_id):
                     enzyme = model.enzymes.get_by_id(enzyme_id)
-
-                    #check if catalytic event already exists:
-                    catal_event_id = f'CE_{rxn_id}'
-
-                    if catal_event_id in model.catalytic_events:
-                        ce = model.catalytic_events.get_by_id(catal_event_id)
-                        enzyme.add_catalytic_event(ce, kcat)
-                    else:
-                        enzyme.create_catalytic_event(rxn_id=rxn_id, kcats=kcat)
-                        model.add_catalytic_events([enzyme.catalytic_events.get_by_id(f'CE_{rxn_id}')])
+                    self._add_reaction_to_enzyme(model, enzyme, rxn_id, kcat)
                 else:
+                    if self.protein2gene != {}:
+                        gene_list = self._get_model_genes_from_enzyme(enzyme_id, model)
+                    else:
+                        gene_list = []
+
                     enzyme = Enzyme(
-                        id = enzyme_id,
-                        rxn2kcat= {rxn_id: kcat},
-                        molmass = molmass
+                        id=enzyme_id,
+                        rxn2kcat={rxn_id: kcat},
+                        rxn2pr={rxn_id: protein_reaction},
+                        molmass=molmass,
+                        genes=gene_list
                     )
+
+                    if self._enzyme_is_enzyme_complex(protein_reaction, enzyme_id):
+                        for pr in protein_reaction:
+                            if len(pr) > 1:
+                                enzyme_complex_id = '_'.join(pr)
+                                if enzyme_complex_id not in model.enzymes:
+                                    enzyme = EnzymeComplex(
+                                        id=enzyme_complex_id,
+                                        rxn2kcat={rxn_id: kcat},
+                                        molmass=0,
+                                        genes=gene_list,
+                                        enzymes=[enzyme]
+                                    )
+                                    model.add_enzymes([enzyme])
+                                    self.constraints += [enzyme]
+                                    self.variables.append(enzyme.enzyme_variable)
+                                else:
+                                    enz_complex = model.enzymes.get_by_id(enzyme_complex_id)
+                                    self._add_reaction_to_enzyme(model, enz_complex, rxn_id, kcat)
+                                    enz_complex.reactions.append(rxn_id)
+                                    enz_complex.add_enzymes([enzyme])
+                    else:
+                        model.add_enzymes([enzyme])
+
+                        self.constraints += [enzyme]
+                        self.variables.append(enzyme.enzyme_variable)
                     #in the add enzymes function the enzyme with corresponding catalytic events will be added to the model
                     #and connected to the reaction and total protein constraint
 
-                    model.add_enzymes([enzyme])
 
                     #adding to the enzyme sector object for easy removal
-                    self.constraints += [enzyme]
                     model.tpc += 1
-                self.variables.append(enzyme.enzyme_variable)
         return model
 
-    def check_kcat_values(self, model, reaction, kcat):
+    def check_kcat_values(self, model, reaction, enzyme_dict):
         """
         Function to check if the kcat values provided for an enzyme are consistent with the direction of the reaction
         Parameters
@@ -176,7 +209,7 @@ def check_kcat_values(self, model, reaction, kcat):
             model to which the kcat values should be added
         reaction: cobra.Reaction
             reaction which is catalyzed with the enzyme related to the kcats
-        kcat: dict
+        enzyme_dict: dict
             a dict with the turnover values for the forward and/or backward reaction for different enzymes [/s]
             {'E1':{'f': forward kcat, 'b': backward kcat}}
 
@@ -193,11 +226,14 @@ def check_kcat_values(self, model, reaction, kcat):
 
          # check consistency between provided kcat values and reaction direction
         directions = []
-        for val in kcat.values():# get all directions from the kcat dict
-            directions += list(val.keys())
         kcats = []
-        for val in kcat.values():# get all directions from the kcat dict
-            kcats += list(val.values())
+        for enzyme_info in enzyme_dict.values():
+            # get all directions from the kcat dict
+            for key, value in enzyme_info.items():
+                if key == 'f' or key =='b':
+                    directions += [key]
+                    kcats += [value]
+
 
         rxn_dir = 'consistent'
         if reaction.lower_bound >= 0 and reaction.upper_bound > 0:
@@ -227,6 +263,52 @@ def check_kcat_values(self, model, reaction, kcat):
 
         return rxn_dir == 'consistent'
 
+    def _add_reaction_to_enzyme(self, model, enzyme:Union[Enzyme, EnzymeComplex], rxn_id:str, kcat:dict):
+        # check if catalytic event already exists:
+        catal_event_id = f'CE_{rxn_id}'
+        if catal_event_id in enzyme.catalytic_events:
+            return
+
+        if catal_event_id in model.catalytic_events:
+            ce = model.catalytic_events.get_by_id(catal_event_id)
+            enzyme.add_catalytic_event(ce, kcat)
+        else:
+            enzyme.create_catalytic_event(rxn_id=rxn_id, kcats=kcat)
+            model.add_catalytic_events([enzyme.catalytic_events.get_by_id(f'CE_{rxn_id}')])
+
+    def _enzyme_is_enzyme_complex(self, protein_reaction, enzyme_id):
+        return any([((enzyme_id in pr) & (len(pr) > 1)) for pr in protein_reaction])
+
+    def _get_model_genes_from_enzyme(self, enzyme_id: str, model: Model) -> list:
+        """
+           Retrieves genes associated with the specified enzyme from the model.
+
+           Args:
+               enzyme_id (str): The identifier of the enzyme.
+               model (Model): The model containing gene information.
+
+           Returns:
+               list: A nested list of gene objects associated with the enzyme.
+           """
+        if enzyme_id not in self.protein2gene.keys(): return []
+        gene_id_list = self.protein2gene[enzyme_id]
+        gene_list = []
+        for genes_or in gene_id_list:
+            genes_and_list = []
+            for gene_and in genes_or:
+                #check if there is an and relation (then gene and should be a list and not a string)
+                if isinstance(gene_and, list):
+                    for gene in gene_and:
+                        if gene not in model.genes:
+                            model.genes.append(Gene(gene))
+                else:
+                    gene = gene_and
+                    if gene not in model.genes:
+                        model.genes.append(Gene(gene))
+                genes_and_list.append(model.genes.get_by_id(gene))
+            gene_list.append(genes_and_list)
+        return gene_list
+
 class TransEnzymeSector(EnzymeSector):
     DEFAULT_MOL_MASS = 4.0590394e05  # default E. coli ribosome molar mass [g/mol]
     BIOMASS_RXNID = Config.BIOMASS_REACTION
@@ -289,6 +371,7 @@ def __init__(self, ups_0, ups_mu ,id_list, mol_mass=None, configuration = Config
             self.ups_mu = ups_mu
         # *1000 to convert units from g/g_cdw to mg/g_cdw
         self.ups_0_coeff = self.ups_0[0]*1e3
+        self.id = 'UnusedEnzymeSector'
 
         self.slope = self.ups_mu*1e3
         self.intercept = self.ups_0_coeff
diff --git a/src/PAModelpy/PAModel.py b/src/PAModelpy/PAModel.py
index b01c28a..3d19372 100644
--- a/src/PAModelpy/PAModel.py
+++ b/src/PAModelpy/PAModel.py
@@ -6,7 +6,7 @@
 #type checking
 from optlang.symbolics import Zero
 from optlang.interface import Objective
-from typing import List, Optional, Union, Dict, Iterable
+from typing import List, Optional, Union, Dict, Iterable, Tuple
 #other
 from functools import partial
 import warnings
@@ -18,7 +18,7 @@
 from .EnzymeSectors import ActiveEnzymeSector, TransEnzymeSector, UnusedEnzymeSector, CustomSector, Sector
 from .CatalyticEvent import CatalyticEvent
 from .Constraints import Constraint
-from .Enzyme import Enzyme
+from .Enzyme import Enzyme, EnzymeComplex
 from .configuration import Config
 
 
@@ -95,7 +95,7 @@ def __init__(self, id_or_model: Union[str, "Model", None] = None,
                  active_sector: Optional[ActiveEnzymeSector]=None,
                  translational_sector: Optional[TransEnzymeSector]=None,
                  unused_sector: Optional[UnusedEnzymeSector]=None,
-                 custom_sectors: Union[List, CustomSector] =None,
+                 custom_sectors: Union[List, CustomSector] =[None],
                  configuration = Config):
         """Constants"""
         self.TOTAL_PROTEIN_CONSTRAINT_ID = configuration.TOTAL_PROTEIN_CONSTRAINT_ID
@@ -132,7 +132,8 @@ def __init__(self, id_or_model: Union[str, "Model", None] = None,
         if sensitivity: #perform sensitivity analysis when the model is run
             self.add_lb_ub_constraints()
 
-        for sector in [active_sector, translational_sector, unused_sector, custom_sectors]:
+        sectors_to_add = [active_sector, translational_sector, unused_sector] + custom_sectors
+        for sector in [sector for sector in sectors_to_add if sector is not None]:
             if sector is not None:
                 self.add_sectors([sector])
         print(f'Done with setting up the proteome allocation model {self.id}\n')
@@ -151,7 +152,7 @@ def translational_enzymes(self):
         return self.sectors.get_by_id('TranslationalEnzymeSector')
 
     @translational_enzymes.setter
-    def tranlational_enzymes(self, slope:float, intercept:float, lin_rxn_id: str = 'BIOMASS_Ec_iML1515_WT_75p37M'):
+    def translational_enzymes(self, slope:float, intercept:float, lin_rxn_id: str = 'BIOMASS_Ec_iML1515_WT_75p37M'):
         #input is in g/gDW
         self.change_sector_parameters(self.translational_enzymes, slope, intercept, lin_rxn_id)
 
@@ -195,6 +196,8 @@ def existing_filter(enz: Enzyme) -> bool:
                 False if enzyme exists, True if it doesn't.
                 If the enzyme exists, will log a warning.
             """
+            if isinstance(enz, EnzymeComplex):
+                return False
             if enz.id in self.enzymes:
                 warnings.warn(f"Ignoring enzyme '{enz.id}' since it already exists.")
                 return False
@@ -273,6 +276,9 @@ def parameter_filter(enz: Enzyme) -> bool:
         if not hasattr(enzyme_list, "__iter__"):
             enzyme_list = [enzyme_list]
         for enz in enzyme_list:
+            # if the enzyme is an enzyme complex, pass it on to the specialized add enzyme complex method
+            if isinstance(enz, EnzymeComplex):
+                self.add_enzyme_complex(enz)
             if not isinstance(enz, Enzyme):
                 raise AttributeError('The input was neither an Iterable nor an Enzyme. Please provide only (lists of) Enzyme objects.')
 
@@ -299,100 +305,120 @@ def parameter_filter(enz: Enzyme) -> bool:
 
         # populate solver
         for enzyme in pruned:
-            #add constraint related to the enzyme
-            self.add_enzyme_constraints([f'EC_{enzyme.id}_f', f'EC_{enzyme.id}_b'])
-            #store the constraints connected to the model in the enzyme object
-            new_constraints = {}
-            new_constraints[f'EC_{enzyme.id}_f'] = self.enzyme_constraints[f'EC_{enzyme.id}_f']
-            new_constraints[f'EC_{enzyme.id}_b'] = self.enzyme_constraints[f'EC_{enzyme.id}_b']
-            enzyme._constraints = new_constraints
-
-            #add enzyme variable
-            enzyme_variable = enzyme.enzyme_variable
-            if enzyme_variable.id not in self.variables:
-                    # make enzyme variable aware of the model.
-                    # The enzyme variable will be added as variable to the model in the model.setter magic function
-                    enzyme_variable.model = self
-                    # self.enzyme_variables.append(enzyme_variable)
-                    if context:
-                        context(partial(setattr, enzyme_variable, "_model", None))
-
-            # connect the enzyme to total protein constraint
-            enzyme_variable = self.enzyme_variables.get_by_id(enzyme.id)
-            self.constraints[self.TOTAL_PROTEIN_CONSTRAINT_ID].set_linear_coefficients(
-                { #*1e-6 to make calculations more accurate (solver accuracy)
+            self._add_enzyme(enzyme)
+
+    def add_enzyme_complex(self, enzyme_complex: EnzymeComplex):
+        if enzyme_complex not in self.enzymes:
+            enzyme_complex.model = self
+            self.enzymes.append(enzyme_complex)
+            self._add_enzyme(enzyme_complex)
+        for enzyme in enzyme_complex.enzymes:
+            if enzyme not in self.enzymes:
+                self.enzymes.append(enzyme)
+
+    def _add_enzyme(self, enzyme: Union[Enzyme, EnzymeComplex]):
+        # add context manager
+        context = get_context(self)
+        # add constraint related to the enzyme
+        self.add_enzyme_constraints([f'EC_{enzyme.id}_f', f'EC_{enzyme.id}_b'])
+        # store the constraints connected to the model in the enzyme object
+        new_constraints = {}
+        new_constraints[f'EC_{enzyme.id}_f'] = self.enzyme_constraints[f'EC_{enzyme.id}_f']
+        new_constraints[f'EC_{enzyme.id}_b'] = self.enzyme_constraints[f'EC_{enzyme.id}_b']
+        enzyme._constraints = new_constraints
+
+        # add enzyme variable
+        enzyme_variable = enzyme.enzyme_variable
+        if enzyme_variable.id not in self.variables:
+            # make enzyme variable aware of the model.
+            # The enzyme variable will be added as variable to the model in the model.setter magic function
+            enzyme_variable.model = self
+            # self.enzyme_variables.append(enzyme_variable)
+            if context:
+                context(partial(setattr, enzyme_variable, "_model", None))
+
+        # connect the enzyme to total protein constraint
+        enzyme_variable = self.enzyme_variables.get_by_id(enzyme.id)
+        self.constraints[self.TOTAL_PROTEIN_CONSTRAINT_ID].set_linear_coefficients(
+            {  # *1e-6 to make calculations more accurate (solver accuracy)
                 enzyme_variable.forward_variable: enzyme.molmass * 1e-6,
                 enzyme_variable.reverse_variable: enzyme.molmass * 1e-6,
-                }
-            )
-            self.tpc+=1
-
-            # add enzyme constraints to the model and link to enzyme and reaction variables
-            #get the enzyme variable
-            enzyme_var_model = self.enzyme_variables.get_by_id(enzyme.id)
-            #connect enzyme variable to its upper and lower bound
-            if enzyme.upper_bound > self.p_tot*1e3:
-                ub = enzyme.upper_bound
+            }
+        )
+        self.tpc += 1
+
+        # add enzyme constraints to the model and link to enzyme and reaction variables
+        # get the enzyme variable
+        enzyme_var_model = self.enzyme_variables.get_by_id(enzyme.id)
+        # connect enzyme variable to its upper and lower bound
+        if enzyme.upper_bound > self.p_tot * 1e3:
+            ub = enzyme.upper_bound
+        else:
+            ub = self.p_tot * 1e3
+        self, enzyme = self.make_enzyme_min_max_constraint(self, enzyme, lower_bound=enzyme.lower_bound, upper_bound=ub)
+
+        # add the enzyme to the interface between the enzyme constraints, enzyme variables and the reaction
+        for catalytic_event in enzyme.catalytic_events:
+            if catalytic_event not in self.catalytic_events:
+                # make catalytic event aware of the model.
+                # The catalytic event will be added configured in the model.setter magic function
+                catalytic_event.model = self
+                self.catalytic_events += [catalytic_event]
+                if context:
+                    context(partial(setattr, catalytic_event, "_model", None))
             else:
-                ub = self.p_tot*1e3
-            self, enzyme = self.make_enzyme_min_max_constraint(self, enzyme, lower_bound=enzyme.lower_bound, upper_bound=ub)
-
-            # add the enzyme to the interface between the enzyme constraints, enzyme variables and the reaction
-            for catalytic_event in enzyme.catalytic_events:
-                if catalytic_event not in self.catalytic_events:
-                    # make catalytic event aware of the model.
-                    # The catalytic event will be added configured in the model.setter magic function
-                    catalytic_event.model = self
-                    self.catalytic_events += [catalytic_event]
-                    if context:
-                        context(partial(setattr, catalytic_event, "_model", None))
-                else:
-                    # add enzyme to the catalytic event
-                    catalytic_event_model = self.catalytic_events.get_by_id(catalytic_event.id)
-                    #remove block on the reaction of no enzyme constraint if this is present
-                    ce_constraints = catalytic_event_model.constraints.copy()
-                    for name, constraint in ce_constraints.items():
-                        if 'no_enzyme' in name:
-                            del catalytic_event_model.constraints[name]
-                            self.remove_cons_vars([constraint])
-
+                # add enzyme to the catalytic event
+                catalytic_event_model = self.catalytic_events.get_by_id(catalytic_event.id)
+                # remove block on the reaction of no enzyme constraint if this is present
+                ce_constraints = catalytic_event_model.constraints.copy()
+                for name, constraint in ce_constraints.items():
+                    if 'no_enzyme' in name:
+                        del catalytic_event_model.constraints[name]
+                        self.remove_cons_vars([constraint])
+
+                if enzyme not in catalytic_event_model.enzymes:
                     catalytic_event_model.add_enzymes({enzyme: enzyme.rxn2kcat[catalytic_event.rxn_id]})
 
-                    # replace the catalytic event with the already existing event from the model
-                    enzyme.catalytic_events._replace_on_id(catalytic_event_model)
-
-            #connect the enzyme to each of the reactions it is associated with
-            for rxn, kcatdict in enzyme.rxn2kcat.items():
-                # link enzymatic variable to reaction via the turnover number kcat
-                for direction, kcat in kcatdict.items():
-                    # check direction
-                    if direction != 'f' and direction != 'b':
-                        warnings.warn(['Invalid kcat direction encountered for ', catalytic_event.id, direction])
-                        continue
-
-                    # create enzyme constraint for the reaction if not existent already
-                    constraint_id = 'EC_' + enzyme.id + '_' + direction
-                    if constraint_id not in self.enzyme_constraints.keys():
-                        self.add_enzyme_constraints([constraint_id])
-
-                    # kcat is used as a coefficient for the enzyme concentration
-                    # Get existent forward/reverse variables
-                    coeff = kcat*3600*1e-6 #3600 to convert to /h to /s *1e-6 to make calculations more accurate
-                    if direction == 'f':
-                        self.constraints[constraint_id].set_linear_coefficients(
-                            {enzyme_var_model.forward_variable: -1,
-                              self.reactions.get_by_id(rxn).forward_variable: 1/coeff
-                            })
-                    elif direction == 'b':
-                        self.constraints[constraint_id].set_linear_coefficients(
-                            {enzyme_var_model.reverse_variable: -1,
-                             self.reactions.get_by_id(rxn).reverse_variable: 1/coeff
-                            })
-
-                    # make reaction-enzyme interface and the enzyme variable aware of its participation in the constraint
-                    catalytic_event_model.constraints[enzyme.id] = self.constraints[constraint_id]
-                    enzyme_var_model.constraints[constraint_id] = self.constraints[constraint_id]
-                    self.solver.update()
+                # replace the catalytic event with the already existing event from the model
+                enzyme.catalytic_events._replace_on_id(catalytic_event_model)
+
+        # connect the enzyme to each of the reactions it is associated with
+        for rxn, kcatdict in enzyme.rxn2kcat.items():
+            # link enzymatic variable to reaction via the turnover number kcat
+            for direction, kcat in kcatdict.items():
+                # check direction
+                if direction != 'f' and direction != 'b':
+                    warnings.warn(['Invalid kcat direction encountered for ', catalytic_event.id, direction])
+                    continue
+
+                # create enzyme constraint for the reaction if not existent already
+                constraint_id = 'EC_' + enzyme.id + '_' + direction
+                if constraint_id not in self.enzyme_constraints.keys():
+                    self.add_enzyme_constraints([constraint_id])
+
+                # kcat is used as a coefficient for the enzyme concentration
+                # Get existent forward/reverse variables
+                self._change_kcat_in_enzyme_constraint(rxn, enzyme.id,
+                                                       direction, kcat)
+                if direction == 'f':
+                    self.constraints[constraint_id].set_linear_coefficients(
+                        {enzyme_var_model.forward_variable: -1,
+                         })
+                elif direction == 'b':
+                    self.constraints[constraint_id].set_linear_coefficients(
+                        {enzyme_var_model.reverse_variable: -1,
+                         })
+
+                # make reaction-enzyme interface and the enzyme variable aware of its participation in the constraint
+                catalytic_event_model.constraints[enzyme.id] = self.constraints[constraint_id]
+                enzyme_var_model.constraints[constraint_id] = self.constraints[constraint_id]
+                self.solver.update()
+
+        # check if all genes are in the model
+        for genes_or in enzyme.genes:
+            for gene_and in genes_or:
+                if not gene_and in self.genes:
+                    self.genes.append(gene_and)
 
     def add_sectors(self, sectors: List = None):
         """
@@ -736,17 +762,17 @@ def make_lb_ub_constraint(m: Optional[Model], rxn: Reaction, lower_bound: float,
     @staticmethod
     def make_enzyme_min_max_constraint(m: Optional[Model], enz: Enzyme, lower_bound: float, upper_bound:float):
         """
-        Adding variables and constraints for the lower and upperbounds of a reaction to a model.
+        Adding variables and constraints for the lower and upperbounds of an Enzyme to a model.
         When solving the model, shadow prices for the lower and upperbounds will be calculated.
         This allows for the calculation of sensitivity coefficients. The constraints are formulated as follows:
-        R_ub : E <= Emax
-        R_lb : -E <= -Emin
+        enz_max : E <= Emax
+        enz_min : -E <= -Emin
 
         Parameters
         ----------
         m: cobra.Model or PAModelpy.PAModel
             model to which the upper and lowerbound constraints and variables should be added
-        rxn: PAModelpy.Enzyme
+        enz: PAModelpy.Enzyme
             Enzyme for which minimal and maximal concentration constraints should be generated
         lower_bound: float
             Value of the lowerbound
@@ -835,19 +861,16 @@ def calculate_csc(self, obj_value, mu, mu_ub, mu_lb, mu_ec_f, mu_ec_b):
 
         Definition: constraint_UB*shadowprice/obj_value.
 
-        :param obj_value: Float
-        :param mu: DataFrame
-            Shadowprices for all constraints
-        :param mu_ub: DataFrame
-            Shadowprices for the reaction UB constraints
-        :param mu_lb: DataFrame
-            Shadowprices for the reaction LB constraints
-        :param mu_ec_f: DataFrame
-            Shadowprices for the constraint related to an enzymatic catalysis of the forward reaction
-        :param mu_ec_b: DataFrame
-            Shadowprices for the constraint related to an enzymatic catalysis of the backward reaction
-
         Results will be saved in the self.capacity_sensitivity_coefficients attribute as a dataframe
+
+        Args:
+            obj_value: Float: optimal objective value, commonly maximal growth rate under specific conditions
+            mu: DataFrame: shadowprices for all constraints
+            mu_ub: DataFrame: Shadowprices for the reaction UB constraints
+            mu_lb: DataFrame: Shadowprices for the reaction LB constraints
+            mu_ec_f: DataFrame: Shadowprices for the constraint related to an enzymatic catalysis of the forward reaction
+            mu_ec_b: DataFrame: Shadowprices for the constraint related to an enzymatic catalysis of the backward reaction
+
         """
         self.capacity_sensitivity_coefficients = pd.DataFrame(columns=['rxn_id', 'enzyme_id', 'constraint', 'coefficient'])
         # add capacity sensitivity coefficients for sectors if they are there
@@ -865,17 +888,14 @@ def calculate_csc(self, obj_value, mu, mu_ub, mu_lb, mu_ec_f, mu_ec_b):
             #treat sectors separately if there is not a total protein constraint
         else:
             for sector in self.sectors:
-                try:
-                    constraint = 'sector'
-                    rxn_id = 'R_' + sector.id
-                    enzyme_id = sector.id
-                    ca_coefficient = self.constraints[enzyme_id].ub * mu[mu['rxn_id'] == sector.id]['shadow_prices'].iloc[0] / obj_value
-
-                    new_row = [rxn_id, enzyme_id, constraint, ca_coefficient]
-                    # add new_row to dataframe
-                    self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_row
-                except:
-                    continue
+                constraint = 'sector'
+                rxn_id = 'R_' + sector.id
+                enzyme_id = sector.id
+                ca_coefficient = self.constraints[enzyme_id].ub * mu[mu['rxn_id'] == sector.id]['shadow_prices'].iloc[0] / obj_value
+
+                new_row = [rxn_id, enzyme_id, constraint, ca_coefficient]
+                # add new_row to dataframe
+                self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_row
 
         for rxn in self.reactions:
             # LB
@@ -895,24 +915,57 @@ def calculate_csc(self, obj_value, mu, mu_ub, mu_lb, mu_ec_f, mu_ec_b):
             self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_row_LB
 
         for enzyme in self.enzymes:
-            # get reactions associated with this enzyme
-            reactions = ','.join(self.get_reactions_with_enzyme_id(enzyme.id))
-            # get the right row from the shadow price dataframes
-            # min constraints can be skipped as they always equal to 0 (the enzymes cannot have a concentration lower than 0)
-            mu_ec_max_row = mu_ec_f[mu_ec_f['index'] == f'{enzyme.id}_max']
-            mu_ec_min_row = mu_ec_b[mu_ec_b['index'] == f'{enzyme.id}_min']
-
-            # max Enzyme constraint
-            ca_coefficient_EC_max= self.constraints[f'{enzyme.id}_max'].ub * mu_ec_max_row['shadow_prices'].iloc[0] / obj_value
-            new_enzyme_row_EC_max =[reactions, enzyme.id, 'enzyme_max', ca_coefficient_EC_max]
-            # add new_row to dataframe
-            self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_enzyme_row_EC_max
+            self.calculate_enzyme_csc(enzyme, mu_ec_f, mu_ec_b, obj_value)
 
-            # min Enzyme constraint
-            ca_coefficient_EC_min = self.constraints[f'{enzyme.id}_min'].ub * mu_ec_min_row['shadow_prices'].iloc[0] / obj_value
-            new_enzyme_row_EC_min =[reactions, enzyme.id, 'enzyme_min', ca_coefficient_EC_min]
-            # add new_row to dataframe
-            self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_enzyme_row_EC_min
+    def calculate_csc_for_molecule(self, molecule: Union[Enzyme],
+                                   mu_min:pd.DataFrame, mu_max:pd.DataFrame, obj_value:float,
+                                   constraint_type:str, associated_reactions:str):
+        """
+        Calculate the capacity sensitivity coefficients (CSCs) for constraints related to a biomolecule,
+        such as enzymes. These coefficients reflect the effect of infitesmal changes in the constraint bounds
+        on the objective function.
+
+        The coefficients and associated reactions will be saved in the capacity_sensitivity_coefficients dataframe.
+
+        Args:
+           enzyme:Enzyme: enzyme object to calculate CSC for
+           mu_min: DataFrame: Shadowprices for the constraint related to a lower bound/minimum
+           mu_max: DataFrame: Shadowprices for the constraint related to an upper bound/maximum
+           obj_value: float: optimal objective value, commonly maximal growth rate under specific conditions
+        """
+        # get the right row from the shadow price dataframes
+        mu_max_row = mu_max[mu_max['index'] == f'{molecule.id}_max']
+        mu_min_row = mu_min[mu_min['index'] == f'{molecule.id}_min']
+
+        # Calculate sensitivity coefficients for maximum constraint
+        if f'{molecule.id}_max' in self.constraints.keys():
+            ca_coefficient_max = self.constraints[f'{molecule.id}_max'].ub * mu_max_row['shadow_prices'].iloc[0] / obj_value
+            new_row_max = [associated_reactions, molecule.id, f'{constraint_type}_max', ca_coefficient_max]
+            self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_row_max
+
+        # Calculate sensitivity coefficients for minimum constraint
+        if f'{molecule.id}_min' in self.constraints.keys():
+            ca_coefficient_min = self.constraints[f'{molecule.id}_min'].ub * mu_min_row['shadow_prices'].iloc[0] / obj_value
+            new_row_min = [associated_reactions, molecule.id, f'{constraint_type}_min', ca_coefficient_min]
+            self.capacity_sensitivity_coefficients.loc[len(self.capacity_sensitivity_coefficients)] = new_row_min
+
+
+    def calculate_enzyme_csc(self, enzyme:Enzyme, mu_ec_f:pd.DataFrame, mu_ec_b:pd.DataFrame, obj_value:float):
+        """
+        Calculate the capacity sensitivity coefficients (CSCs) for constraints related to enzyme. These coefficients
+        reflect the effect of infitesmal changes in the constraint bounds on the objective function. The coefficients
+        and associated reactions will be saved in the capacity_sensitivity_coefficients dataframe.
+
+        The function makes use of the abstracted function calculate_csc_for_molecule
+
+        Args:
+            enzyme:Enzyme: enzyme object to calculate CSC for
+            mu_ec_f: DataFrame: Shadowprices for the constraint related to an enzymatic catalysis of the forward reaction
+            mu_ec_b: DataFrame: Shadowprices for the constraint related to an enzymatic catalysis of the backward reaction
+            obj_value: float: optimal objective value, commonly maximal growth rate under specific conditions
+        """
+        reactions = ','.join(self.get_reactions_with_enzyme_id(enzyme.id))
+        self.calculate_csc_for_molecule(enzyme, mu_ec_b, mu_ec_f, obj_value, 'enzyme', reactions)
 
     def calculate_esc(self, obj_value, mu_ec_f, mu_ec_b):
         """
@@ -1059,6 +1112,34 @@ def change_reaction_lb(self,rxn_id:str, lower_bound: float = None):
         else:
             self.reactions.get_by_id(rxn_id).lower_bound = lower_bound
 
+    def get_reaction_bounds(self, rxn_id:str) -> Tuple[Union[float, int]]:
+        """
+        Get the reaction bounds. If there should be a sensitivity analysis, the bounds of the upper and lower bound
+        constraints returned
+        Args:
+        rxn_id: str
+            string of reaction id to return
+        """
+        if rxn_id not in self.reactions:
+            warnings.warn(f'Reaction {rxn_id} does not exist in the model. Cannot get the upper- and lowerbound.')
+            return
+        #make sure the order of setting is right to prevent errors
+        lb = self.get_reaction_lb(rxn_id)
+        ub = self.get_reaction_ub(rxn_id)
+        return lb, ub
+
+    def get_reaction_ub(self, rxn_id:str) -> Union[int, float]:
+        if self.sensitivity:
+            return self.constraints[rxn_id + '_ub'].ub
+        else:
+            return self.reactions.get_by_id(rxn_id).upper_bound
+
+    def get_reaction_lb(self,rxn_id:str):
+        if self.sensitivity:
+            return -self.constraints[rxn_id + '_lb'].ub
+        else:
+            return self.reactions.get_by_id(rxn_id).lower_bound
+
     def change_enzyme_bounds(self, enzyme_id: str, lower_bound: float = None, upper_bound: float = None):
         """
                 Change the enzyme bounds. If the model should be primed for performing a sensitivity analysis,
@@ -1097,7 +1178,7 @@ def change_enzyme_min(self, enzyme_id: str, lower_bound: float = None):
         else:
             self.enzyme_variables.get_by_id(enzyme_id).lower_bound = lower_bound
 
-    def get_enzymes_with_reaction_id(self, rxn_id:str):
+    def get_enzymes_with_reaction_id(self, rxn_id:str) -> DictList:
         """
         Returns Enzyme objects associated with the reaction id through CatalyticEvent objects
         :param rxn_id: str
@@ -1107,6 +1188,9 @@ def get_enzymes_with_reaction_id(self, rxn_id:str):
         -------
             DictList of Enzyme objects associated with the reaction
         """
+        if rxn_id not in self.reactions:
+            warnings.warn(f'Reaction {rxn_id} is not in the model')
+            return DictList()
         catalytic_event_id = 'CE_' + rxn_id
         catalytic_event = self.catalytic_events.get_by_id(catalytic_event_id)
         enzymes = catalytic_event.enzymes
@@ -1143,6 +1227,26 @@ def change_kcat_value(self, enzyme_id:str, kcats:dict):
         else:
             warnings.warn(f'The enzyme {enzyme_id} does not exist in the model. The kcat can thus not be changed.')
 
+    def _change_kcat_in_enzyme_constraint(self, rxn:Union[str, cobra.Reaction], enzyme_id: str,
+                                                direction: str, kcat: float):
+        constraint_id = f'EC_{enzyme_id}_{direction}'
+        if isinstance(rxn, str):
+            rxn = self.reactions.get_by_id(rxn)
+        # change kcat value in the constraint
+        if kcat == 0:
+            coeff = 0
+        else:
+            coeff = 1 / (kcat * 3600 * 1e-6) #3600 to convert to /h to /s *1e-6 to make calculations more accurate
+        if direction == 'f':
+            self.constraints[constraint_id].set_linear_coefficients({
+                rxn.forward_variable: coeff
+            })
+        elif direction == 'b':
+            self.constraints[constraint_id].set_linear_coefficients({
+                rxn.reverse_variable: coeff
+            })
+        self.solver.update()
+
     def remove_enzymes(self,
                        enzymes: Union[str, Enzyme, List[Union[str, Enzyme]]]
                        )-> None:
@@ -1193,6 +1297,48 @@ def remove_enzymes(self,
             for group in associated_groups:
                 group.remove_members(enzyme)
 
+    def remove_enzyme_reaction_association(self,
+                                           enzyme: Union[str, Enzyme],
+                                           reaction: Union[str, Reaction],
+                                           )-> None:
+        """Remove an enzyme-reaction association from the model. Adapted from the cobra.core.remove_reactions() function.
+        If the reaction is not catalyzed by any enzyme anymore, the reaction ub will become 0
+
+                Args:
+                enzyme : Enzyme or str
+                    An enzyme, or the enzyme id for which the association should be removed to remove.
+                reaction : Reaction or str
+                    A reaction, or the reaction id for which the association should be removed to remove.
+
+                """
+
+        if isinstance(reaction, str): reaction = self.reactions.get_by_id(reaction)
+
+        try:
+            enzyme = self.enzymes[self.enzymes.index(enzyme)]
+        except ValueError:
+            warnings.warn(f"{enzyme.id} not in {self}")
+
+        #remove enzyme from catalytic event
+        catalytic_event = self.catalytic_events.get_by_id('CE_' + reaction.id)
+
+        # removing catalytic event from the enzymes
+        for enzyme in catalytic_event.enzymes:
+            if catalytic_event in enzyme.catalytic_events:
+                enzyme.remove_catalytic_event(catalytic_event)
+
+        for enzyme_var in catalytic_event.enzyme_variables:
+            if catalytic_event in enzyme_var.catalytic_events:
+                enzyme_var.remove_catalytic_event(catalytic_event)
+
+        #remove reaction from enzyme constraint
+        for constraint in enzyme._constraints.values():
+            constraint.set_linear_coefficients(
+                {reaction.forward_variable:0,
+                 reaction.reverse_variable:0
+                }
+            )
+
     def remove_reactions(
         self,
         reactions: Union[str, Reaction, List[Union[str, Reaction]]],
@@ -1506,7 +1652,8 @@ def pfba(self,
         self.objective.set_linear_coefficients({v: 1.0 for v in variables})
 
         #run pFBA
-        self.optimize()
+        sol = self.optimize()
+        return sol
 
     def reset_objective(self):
         """
diff --git a/tests/unit_tests/test_pamodel/test_pam_generation.py b/tests/unit_tests/test_pamodel/test_pam_generation.py
new file mode 100644
index 0000000..2a0fac1
--- /dev/null
+++ b/tests/unit_tests/test_pamodel/test_pam_generation.py
@@ -0,0 +1,194 @@
+import pytest
+from src.PAModelpy.configuration import Config
+from src.PAModelpy.PAModel import PAModel
+
+
+from Scripts.toy_ec_pam import build_toy_gem, build_active_enzyme_sector, build_translational_protein_sector, build_unused_protein_sector
+from Scripts.pam_generation_new import (set_up_ecolicore_pam, set_up_ecoli_pam, set_up_toy_pam,
+                                        parse_gpr_information_for_protein2genes,
+                                        parse_gpr_information_for_rxn2protein)
+
+
+def test_gpr_information_is_parsed_correctly():
+    # Arrange
+    enzyme_id = '1.1.5.3'
+    gpr_string ='(b0902 and b0903) or (b0902 and b3114) or (b3951 and b3952) or ((b0902 and b0903) and b2579)'
+
+    parsed_gpr_reference = [['b0902','b0903'], ['b0902','b3114'], ['b3951' ,'b3952'],['b0902', 'b0903','b2579']]
+
+    # Act
+    gpr_parsed = parse_gpr_information_for_protein2genes(gpr_string)
+
+    # Assert
+    assert gpr_parsed == parsed_gpr_reference
+
+def test_gpr_information_for_protein_is_correctly_filtered():
+    # Arrange
+    enzyme_id = 'P0AFH2'
+    gpr_string ='((b1245 and b1247 and b1244 and b1329 and b1246) or (b2177 and b2180 and b2178 and b2179) or (b1245 and b1247 and b1246 and b1244 and b1243))'
+    protein2gene = {'P02932':'b0241', 'P77348': 'b1329','P76027': 'b1246','P33913': 'b2177','P77737': 'b1247',
+                    'P0AFU0': 'b2178', 'P33916': 'b2180', 'P0AFH6': 'b1245', 'P23843': 'b1243','P0AFH2': 'b1244',
+                    'P33915': 'b2179'}
+    gene2protein = {v: k for k, v in protein2gene.items()}
+
+    parsed_gr_reference = [['b1245','b1247','b1244','b1329','b1246'],
+                           ['b1245','b1247','b1246','b1244','b1243']]
+
+    parsed_pr_reference = [['P0AFH6', 'P77737', 'P0AFH2', 'P77348', 'P76027'],
+                           ['P0AFH6', 'P77737', 'P76027', 'P0AFH2', 'P23843']]
+
+    # Act
+    gr_parsed, pr_parsed = parse_gpr_information_for_rxn2protein(gpr_string, gene2protein, protein2gene, enzyme_id)
+
+    # Assert
+    assert gr_parsed == parsed_gr_reference
+    assert pr_parsed == parsed_pr_reference
+
+def test_if_enzyme_complex_in_toy_pam_is_parsed_correctly():
+    sut = set_up_toy_pam_with_enzyme_complex(sensitivity=False)
+
+    assert all([enz in sut.enzymes for enz in ['E1', 'E2', 'E10', 'E2_E10']])
+    assert all([const not in sut.constraints.keys() for const in ['EC_E10_f', 'EC_E2_f']])
+    constraint = sut.constraints['EC_E2_E10_f'].get_linear_coefficients([sut.reactions.CE_R2_E2_E10.forward_variable])
+    assert constraint[sut.reactions.CE_R2_E2_E10.forward_variable] > 0
+
+def test_if_isozymes_in_toy_pam_are_parsed_correctly():
+    sut = set_up_toy_pam_with_isozymes(sensitivity=False)
+
+    #check whether catalytic reactions are configured correclty
+    neg_constraint = sut.constraints['CE_R2'].get_linear_coefficients([sut.reactions.R2.reverse_variable,
+                                                                   sut.reactions.CE_R2_E2.forward_variable,
+                                                                   sut.reactions.CE_R2_E10.forward_variable])
+
+    pos_constraint = sut.constraints['CE_R2'].get_linear_coefficients([sut.reactions.R2.forward_variable,
+                                                                   sut.reactions.CE_R2_E2.reverse_variable,
+                                                                   sut.reactions.CE_R2_E10.reverse_variable])
+
+    assert all([enz in sut.enzymes for enz in ['E1', 'E2', 'E10']])
+    assert all([const in sut.constraints.keys() for const in ['EC_E10_f', 'EC_E2_f']])
+    assert all([rxn in sut.reactions for rxn in ['CE_R1_E1', 'CE_R2_E2', 'CE_R2_E10']])
+    assert all([v == -1 for v in neg_constraint.values()])
+    assert all([v == 1 for v in pos_constraint.values()])
+
+def test_if_toy_pam_with_isozymes_has_same_growth_rate_as_without():
+    sut = set_up_toy_pam_with_isozymes(sensitivity=False)
+    toy_pam = set_up_toy_pam()
+
+    for model in [sut, toy_pam]:
+        model.change_reaction_bounds('R1', 0,0.01)
+        model.optimize()
+
+    assert sut.objective.value == pytest.approx(toy_pam.objective.value, abs = 1e-6)
+
+def test_if_toy_pam_with_isozymes_and_enzymecomplex_has_same_growth_rate_as_without():
+    sut = set_up_toy_pam_with_isozymes_and_enzymecomplex(sensitivity=False)
+    toy_pam = set_up_toy_pam()
+
+    for model in [sut, toy_pam]:
+        model.change_reaction_bounds('R1', 0,0.01)
+        model.optimize()
+
+    assert sut.objective.value == pytest.approx(toy_pam.objective.value, abs = 1e-6)
+
+
+def test_if_toy_pam_with_enzyme_comples_has_same_growth_rate_as_without():
+    sut = set_up_toy_pam_with_enzyme_complex(sensitivity=False)
+    toy_pam = set_up_toy_pam()
+
+    for model in [sut, toy_pam]:
+        model.change_reaction_bounds('R1', 0,0.01)
+        model.optimize()
+
+    assert sut.objective.value == pytest.approx(toy_pam.objective.value, abs = 1e-6)
+
+# def test_set_up_ecolicore_pam_works():
+#     sut = set_up_ecolicore_pam()
+#     assert True
+# def test_set_up_ecoli_pam_works():
+#     sut = set_up_ecoli_pam()
+#     assert True
+
+#########################################################################################################################
+# HELPER FUNCTIONS
+##########################################################################################################################
+
+def set_up_toy_pam_with_enzyme_complex(sensitivity =True):
+    config = Config()
+    #setting the configuration for the toy model
+    config.BIOMASS_REACTION = 'R7'
+    config.GLUCOSE_EXCHANGE_RXNID = 'R1'
+    config.CO2_EXHANGE_RXNID = 'R8'
+    config.ACETATE_EXCRETION_RXNID = 'R9'
+
+    Etot = 0.6*1e-3
+    model = build_toy_gem()
+    active_enzyme = build_active_enzyme_sector(config)
+    #add an enzyme associated to enzyme complex to the toy model
+    active_enzyme.rxn2protein['R2']['E2']['protein_reaction_association'] = [['E2', 'E10']]
+    active_enzyme.rxn2protein['R2']['E10']= active_enzyme.rxn2protein['R2']['E2'].copy()
+
+    #build the toy model
+    unused_enzyme = build_unused_protein_sector(config)
+    translation_enzyme = build_translational_protein_sector(config)
+    pamodel = PAModel(model, name='toy model MCA with enzyme constraints', active_sector=active_enzyme,
+                      translational_sector=translation_enzyme,
+                      unused_sector=unused_enzyme, p_tot=Etot, sensitivity=sensitivity)
+    pamodel.objective = 'R7'
+    config.reset()
+    return pamodel
+
+def set_up_toy_pam_with_isozymes(sensitivity =True):
+    config = Config()
+    #setting the configuration for the toy model
+    config.BIOMASS_REACTION = 'R7'
+    config.GLUCOSE_EXCHANGE_RXNID = 'R1'
+    config.CO2_EXHANGE_RXNID = 'R8'
+    config.ACETATE_EXCRETION_RXNID = 'R9'
+
+    Etot = 0.6*1e-3
+    model = build_toy_gem()
+    active_enzyme = build_active_enzyme_sector(config)
+    #add an enzyme associated to isozymes to the toy model
+    active_enzyme.rxn2protein['R2']['E2']['protein_reaction_association'] = [['E2'], ['E10']]
+    active_enzyme.rxn2protein['R2']['E10']= active_enzyme.rxn2protein['R2']['E2'].copy()
+
+    #build the toy model
+    unused_enzyme = build_unused_protein_sector(config)
+    translation_enzyme = build_translational_protein_sector(config)
+    pamodel = PAModel(model, name='toy model MCA with enzyme constraints', active_sector=active_enzyme,
+                      translational_sector=translation_enzyme,
+                      unused_sector=unused_enzyme, p_tot=Etot, sensitivity=sensitivity)
+    pamodel.objective = 'R7'
+    config.reset()
+    return pamodel
+
+def set_up_toy_pam_with_isozymes_and_enzymecomplex(sensitivity =True):
+    config = Config()
+    #setting the configuration for the toy model
+    config.BIOMASS_REACTION = 'R7'
+    config.GLUCOSE_EXCHANGE_RXNID = 'R1'
+    config.CO2_EXHANGE_RXNID = 'R8'
+    config.ACETATE_EXCRETION_RXNID = 'R9'
+
+    Etot = 0.6*1e-3
+    model = build_toy_gem()
+    active_enzyme = build_active_enzyme_sector(config)
+    #add an enzyme associated to isozymes to the toy model
+    active_enzyme.rxn2protein['R2']['E2']['protein_reaction_association'] = [['E2'], ['E10']]
+    active_enzyme.rxn2protein['R2']['E10']= active_enzyme.rxn2protein['R2']['E2'].copy()
+
+    #add an enzyme associated to enzyme complex to the toy model
+    active_enzyme.rxn2protein['R3']['E3']['protein_reaction_association'] = [['E3','E10', 'E11']]
+    active_enzyme.rxn2protein['R3']['E10']= active_enzyme.rxn2protein['R3']['E3'].copy()
+    active_enzyme.rxn2protein['R3']['E11']= active_enzyme.rxn2protein['R3']['E3'].copy()
+
+
+    #build the toy model
+    unused_enzyme = build_unused_protein_sector(config)
+    translation_enzyme = build_translational_protein_sector(config)
+    pamodel = PAModel(model, name='toy model MCA with enzyme constraints', active_sector=active_enzyme,
+                      translational_sector=translation_enzyme,
+                      unused_sector=unused_enzyme, p_tot=Etot, sensitivity=sensitivity)
+    pamodel.objective = 'R7'
+    config.reset()
+    return pamodel