Add total_rate option to rules for NFsim compatibility (#559)

* Add total_rate option to rules for NFSim compatibility
* Test for total_rate functionality in NFsim

pysb/core.py

@@ -1432,6 +1432,13 @@ class Rule(Component):
         co-transport anything connected to that Monomer by a path in the same
         compartment. If False (default), connected Monomers will remain where
         they were.
+    total_rate: bool, optional
+        If True, the rate is considered to be macroscopic and is not
+        multiplied by the number of reactant molecules during simulation.
+        If False (default), the rate is multiplied by number of reactant
+        molecules.
+        Keyword is used by BioNetGen only for simulations using NFsim.
+        Keyword is ignored by generate_network command of BioNetGen.
 
     Attributes
     ----------
@@ -1443,7 +1450,7 @@ class Rule(Component):
 
     def __init__(self, name, rule_expression, rate_forward, rate_reverse=None,
                  delete_molecules=False, move_connected=False,
-                 _export=True):
+                 _export=True, total_rate=False):
         if not isinstance(rule_expression, RuleExpression):
             raise Exception("rule_expression is not a RuleExpression object")
         validate_expr(rate_forward, "forward rate")
@@ -1460,6 +1467,7 @@ def __init__(self, name, rule_expression, rate_forward, rate_reverse=None,
         self.rate_reverse = rate_reverse
         self.delete_molecules = delete_molecules
         self.move_connected = move_connected
+        self.total_rate = total_rate
         # TODO: ensure all numbered sites are referenced exactly twice within each of reactants and products
 
         # Check synthesis products are concrete

pysb/generator/bng.py

@@ -102,6 +102,8 @@ def generate_reaction_rules(self):
                 self.__content += ' DeleteMolecules'
             if r.move_connected:
                 self.__content += ' MoveConnected'
+            if r.total_rate:
+                self.__content += ' TotalRate'
             self.__content += "\n"
         self.__content += "end reaction rules\n\n"
 

pysb/tests/test_bng.py

@@ -156,6 +156,55 @@ def test_nfsim():
         assert len(res) == 101
 
 
+@with_model
+def test_nfsim_total_rate():
+    Monomer('A', ['a'])
+    Monomer('B', ['b'])
+
+    Parameter('ksynthA', 100)
+    Parameter('ksynthB', 100)
+    Parameter('kbindAB', 10)
+
+    Parameter('A_init', 20)
+    Parameter('B_init', 30)
+
+    Initial(A(a=None), A_init)
+    Initial(B(b=None), B_init)
+
+    Observable("A_free", A(a=None))
+    Observable("B_free", B(b=None))
+    Observable("AB_complex", A(a=1) % B(b=1))
+
+    Rule('A_synth', None >> A(a=None), ksynthA)
+    Rule('B_synth', None >> B(b=None), ksynthB)
+    Rule('AB_bind', A(a=None) + B(b=None) >> A(a=1) % B(b=1), kbindAB)
+
+    with BngFileInterface(model) as bng:
+        bng.action('simulate', method='nf', t_end=1000, n_steps=100)
+        bng.execute()
+        res = bng.read_simulation_results()
+        assert res.dtype.names == ('time', 'A_free', 'B_free', 'AB_complex')
+        no_total_rate_ab_complex = res['AB_complex'][-1]
+
+    # Set total rate of last AB_bind reaction to be constant and independent
+    # of A and B concentration. In this case, the number of AB complex molecules
+    # should be less than the non total rate case.
+    model.rules[-1].total_rate = True
+
+    # check that generate network does not fail when total rate is set to be true
+    # generate_network just ignores this setting
+    ok_(generate_network(model))
+
+    with BngFileInterface(model) as bng:
+        bng.action('simulate', method='nf', t_end=1000, n_steps=100)
+        bng.execute()
+        res = bng.read_simulation_results()
+        assert res.dtype.names == ('time', 'A_free', 'B_free', 'AB_complex')
+        total_rate_ab_complex = res['AB_complex'][-1]
+
+    assert no_total_rate_ab_complex > total_rate_ab_complex
+
+
 @with_model
 def test_unicode_strs():
     Monomer(u'A', [u'b'], {u'b':[u'y', u'n']})