diff --git a/xbanalysis/model/fbc_gene_product.py b/xbanalysis/model/fbc_gene_product.py
new file mode 100644
index 0000000000000000000000000000000000000000..efd547704b0f24e93b972cd059221cafac9b9614
--- /dev/null
+++ b/xbanalysis/model/fbc_gene_product.py
@@ -0,0 +1,12 @@
+"""Implementation of FbcGeneProduct class.
+
+Peter Schubert, HHU Duesseldorf, July 2022
+"""
+
+
+class FbcGeneProduct:
+
+    def __init__(self, s_fbc_gene_product):
+        self.id = s_fbc_gene_product.name
+        self.name = s_fbc_gene_product.get('name', '')
+        self.label = s_fbc_gene_product['label']
diff --git a/xbanalysis/model/sbo_term.py b/xbanalysis/model/sbo_term.py
index 8272d072a9fb5c99d282a43ac6a465baa56384ad..52b5ace86069be7e01c5fea198dd4b20a66e8e6e 100644
--- a/xbanalysis/model/sbo_term.py
+++ b/xbanalysis/model/sbo_term.py
@@ -97,7 +97,7 @@ class SboTerm:
 
     @property
     def sbo_term(self):
-        return self.__sbo_term
+        return re.sub('_', ':', self.__sbo_term)
 
     @sbo_term.setter
     def sbo_term(self, sbo_term):
diff --git a/xbanalysis/model/xba_model.py b/xbanalysis/model/xba_model.py
index 3f643b91d4d78fd51cb7b052768412c9161cb04f..29bcbccbd4b531cd1374efb9c81b6753d9f595eb 100644
--- a/xbanalysis/model/xba_model.py
+++ b/xbanalysis/model/xba_model.py
@@ -14,6 +14,7 @@ from .xba_parameter import XbaParameter
 from .xba_reaction import XbaReaction
 from .xba_species import XbaSpecies
 from .fbc_objective import FbcObjective
+from .fbc_gene_product import FbcGeneProduct
 
 
 # TODO: implement isGBAmodel(), isFBAmodel(), isRBAmodel()
@@ -51,6 +52,9 @@ class XbaModel:
         if 'fbcObjectives' in model_dict:
             self.fbc_objectives = {oid: FbcObjective(row)
                                    for oid, row in model_dict['fbcObjectives'].iterrows()}
+        if 'fbcGeneProducts' in model_dict:
+            self.fbc_gene_products = {gp_id: FbcGeneProduct(row)
+                                      for gp_id, row in model_dict['fbcGeneProducts'].iterrows()}
 
     def get_stoic_matrix(self, sids, rids):
         """retrieve stoichiometric sub-matrix [sids x rids].
diff --git a/xbanalysis/model/xba_reaction.py b/xbanalysis/model/xba_reaction.py
index 7bbc2b6a91966d0158da4532242de51ee835c9fd..e483fee3fccab608621b74d3038b5c5433a58f4d 100644
--- a/xbanalysis/model/xba_reaction.py
+++ b/xbanalysis/model/xba_reaction.py
@@ -47,6 +47,9 @@ class XbaReaction:
             if self.reversible is False:
                 self.fbc_lb = max(0, self.fbc_lb)
                 self.fbc_ub = max(0, self.fbc_ub)
+        if 'fbcGeneProdAssoc' in s_reaction:
+            if type(s_reaction['fbcGeneProdAssoc']) == str:
+                self.fbc_gpa = re.sub(r'^assoc=', '', s_reaction['fbcGeneProdAssoc'])
         if 'xmlAnnotation' in s_reaction:
             attrs = sbmlxdf.extract_xml_attrs(s_reaction['xmlAnnotation'], ns=xml_rba_ns, token='kinetics')
             if 'kapp_f_per_s' in attrs:
diff --git a/xbanalysis/model/xba_species.py b/xbanalysis/model/xba_species.py
index a7afc6243ca32b103050d35cd6f63842f9794459..71706e7cd7f43d3da3d8dd9201a765ce19fe7f34 100644
--- a/xbanalysis/model/xba_species.py
+++ b/xbanalysis/model/xba_species.py
@@ -34,7 +34,7 @@ class XbaSpecies:
             if 'process_capacity_per_s' in attrs:
                 self.rba_process_capacity = float(attrs.pop('process_capacity_per_s'))
             if len(attrs) > 0:
-                self.rba_process_costs = attrs
+                self.rba_process_costs = {key: float(val) for key, val in attrs.items()}
 
         self.ps_rid = None
         self.m_reactions = {}
diff --git a/xbanalysis/problems/fba_problem.py b/xbanalysis/problems/fba_problem.py
index ed73a4b96a57819656ddfb09e67d624862bc6279..539a47a580e295352d6523897ae79884404c0538 100644
--- a/xbanalysis/problems/fba_problem.py
+++ b/xbanalysis/problems/fba_problem.py
@@ -143,7 +143,6 @@ class FbaProblem:
         if lp is None:
             return {'success': False, 'message': 'not an FBA model'}
         res = lp.solve()
-        del lp
         res['x'] = self.combine_fwd_bwd(res['x'])
         res['reduced_costs'] = self.combine_fwd_bwd(res['reduced_costs'])
         return res
@@ -176,7 +175,6 @@ class FbaProblem:
             res = lp.solve()
             res['x'] = self.combine_fwd_bwd(res['x'])
             res['reduced_costs'] = self.combine_fwd_bwd(res['reduced_costs'])
-        del lp
         return res
 
     def fva_optimize(self, rids=None, fract_of_optimum=1.0):
@@ -231,7 +229,6 @@ class FbaProblem:
                 else:
                     print(f"FVA max failed for {rid}")
                     flux_ranges[1, idx] = np.nan
-        del lp
 
         res = {'rids': rids, 'min': flux_ranges[0], 'max': flux_ranges[1],
                'obj_min': flux_ranges[2], 'obj_max': flux_ranges[3]}
diff --git a/xbanalysis/problems/rba_problem.py b/xbanalysis/problems/rba_problem.py
index 13433e9a0d679b3c00b08a2dca98fbdfeaeb8284..11c6dd7ecfd7686a1cc95e736fcbbb5d522e5b0b 100644
--- a/xbanalysis/problems/rba_problem.py
+++ b/xbanalysis/problems/rba_problem.py
@@ -24,6 +24,8 @@ class RbaProblem:
 
         :param xba_model:
         """
+        self.xba_model = xba_model
+
         # metabolic reactions and enzyme transitions (no FBA reactions, no protein synthesis, no degradation)
         self.mr_rids = [r.id for r in xba_model.reactions.values() if r.sboterm.is_in_branch('SBO:0000167')
                         and r.sboterm.is_in_branch('SBO:0000179') is False]
@@ -145,7 +147,6 @@ class RbaProblem:
             lp = LpProblem()
             lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
             res = lp.solve(scaling=self.scaling)
-            del lp
             if res['success'] and res['fun'] > 1e-10:
                 gr['lb'] = gr['try']
                 gr['try'] *= 10
@@ -169,13 +170,19 @@ class RbaProblem:
             lp = LpProblem()
             lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
             res = lp.solve(scaling=self.scaling)
-            del lp
             if res['success'] and res['fun'] > 1e-10:
                 gr['lb'] = gr['try']
             else:
                 gr['ub'] = gr['try']
             if (gr['ub'] - gr['lb']) < gr['try'] * 1e-5:
                 break
+
+        # get optimum values for the final solution
+        gr['try'] = gr['lb']
+        constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+        lp = LpProblem()
+        lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+        res = lp.solve(scaling=self.scaling)
         res['fun'] = gr['try']
         return res