Make p_bound function unique and shared across rate laws

tests/eukaryote/test_transcription.py

@@ -267,18 +267,18 @@ def test_methods(self):
         self.assertEqual(model.parameters.get_one(id='total_mitochondrial_genome_binding').units, unit_registry.parse_units('molecule'))
         self.assertEqual(model.parameters.get_one(id='total_mitochondrial_genome_binding').comments, 'Set to genome length divided by 2 bp')
         self.assertEqual(model.parameters.get_one(id='total_mitochondrial_genome_binding').references[0].pages, '394-403')
-        self.assertEqual(model.functions.get_one(id='p_bound_gene1').expression.expression,
+        self.assertEqual(model.functions.get_one(id='p_bound_1').expression.expression,
             '1 / (1 + total_nuclear_genome_binding / (total_complex1_n * 1) * exp(log(K_d_specific_polr / K_d_non_specific_polr)))')
-        self.assertEqual(model.functions.get_one(id='p_bound_gene2').expression.expression,
+        self.assertEqual(model.functions.get_one(id='p_bound_2').expression.expression,
             '1 / (1 + total_mitochondrial_genome_binding / (total_complex3_m * ((1 + activator[m] / (Ka_transcription_initiation_trans2_activator * '
             'Avogadro * volume_m) * f_transcription_initiation_trans2_activator) / (1 + activator[m] / '
             '(Ka_transcription_initiation_trans2_activator * Avogadro * volume_m))) * (1 / (1 + repressor[m] / '
             '(Kr_transcription_initiation_trans2_repressor * Avogadro * volume_m)))) * '
             'exp(log(K_d_specific_polr / K_d_non_specific_polr)))')
         self.assertEqual(model.rate_laws.get_one(id='transcription_initiation_trans1-forward').expression.expression,
-            'p_bound_gene1 * k_specific_binding_complex1 * complex1_bound_non_specific_site[n] * max(min(gene1_binding_site[n] , 1) , 0)')
+            'p_bound_1 * k_specific_binding_complex1 * complex1_bound_non_specific_site[n] * max(min(gene1_binding_site[n] , 1) , 0)')
         self.assertEqual(model.rate_laws.get_one(id='transcription_initiation_trans2-forward').expression.expression,
-            'p_bound_gene2 * k_specific_binding_complex3 * complex3_bound_non_specific_site[m] * max(min(gene2_binding_site[m] , 1) , 0)')
+            'p_bound_2 * k_specific_binding_complex3 * complex3_bound_non_specific_site[m] * max(min(gene2_binding_site[m] , 1) , 0)')
 
         # elongation
         self.assertEqual(model.rate_laws.get_one(id='transcription_elongation_trans2-forward').expression.expression,

wc_model_gen/eukaryote/transcription.py

@@ -536,6 +536,8 @@ def gen_rate_laws(self):
 
         # Generate rate laws for initiation and elongation & termination
         polr_occupancy_width = self.options.get('polr_occupancy_width')
+        p_function_exprs = {}
+        self._gene_p_function_map = {}
         rate_law_no = 0                    
         rnas_kb = cell.species_types.get(__type=wc_kb.eukaryote.TranscriptSpeciesType)
         for rna_kb in rnas_kb:
@@ -619,21 +621,27 @@ def gen_rate_laws(self):
                 Kd_specific_polr.id,
                 Kd_non_specific_polr.id
                 )
-            p_bound_expression, error = wc_lang.FunctionExpression.deserialize(p_bound, {
-                wc_lang.Species: reg_species,
-                wc_lang.Parameter: reg_parameters,            
-                wc_lang.Function: reg_functions,
-                wc_lang.Observable: {polr_obs.id: polr_obs}
-                })
-            assert error is None, str(error)
-
-            p_bound_function = model.functions.create(
-                id='p_bound_{}'.format(rna_kb.gene.id),
-                name='probability of RNAP binding to {}'.format(rna_kb.gene.name),
-                expression=p_bound_expression,
-                references=[ref_model],
-                units=unit_registry.parse_units(''),
-                )
+            if p_bound in p_function_exprs:
+                p_bound_function = p_function_exprs[p_bound]
+                self._gene_p_function_map[rna_kb.id] = p_bound_function
+            else:    
+                p_bound_expression, error = wc_lang.FunctionExpression.deserialize(p_bound, {
+                    wc_lang.Species: reg_species,
+                    wc_lang.Parameter: reg_parameters,            
+                    wc_lang.Function: reg_functions,
+                    wc_lang.Observable: {polr_obs.id: polr_obs}
+                    })
+                assert error is None, str(error)
+
+                p_bound_function = model.functions.create(                    
+                    name='probability of RNAP binding to {}'.format(rna_kb.gene.name),
+                    expression=p_bound_expression,
+                    references=[ref_model],
+                    units=unit_registry.parse_units(''),
+                    )
+                p_bound_function.id = 'p_bound_{}'.format(len(p_function_exprs)+1)
+                p_function_exprs[p_bound] = p_bound_function
+                self._gene_p_function_map[rna_kb.id] = p_bound_function
 
             specific_binding_constant = model.parameters.get_one(
                 id='k_specific_binding_{}'.format(polr_complex_species.species_type.id))
@@ -753,7 +761,7 @@ def calibrate_submodel(self):
             for i in self._polr_pool[rna_pol_pair[rna_kb.id]]:
                 init_reg_species_count[i] = total_polr / no_of_polr_pool            
 
-            p_bound_function = model.functions.get_one(id='p_bound_{}'.format(rna_kb.gene.id))
+            p_bound_function = self._gene_p_function_map[rna_kb.id]
             p_bound_value = p_bound_function.expression._parsed_expression.eval({
                 wc_lang.Species: init_reg_species_count,
                 wc_lang.Compartment: {