Bug in 'bind_table_complex' macro

[lh64]

The example given for the bind_table_complex macro is

bind_table_complex([[C1(y=1) % C2(y=1), C2],
[R1, (1e-4, 1e-1), (2e-4, 2e-1)],
[R1(c1=1) % R2(x=1), (3e-4, 3e-1), None]], 'x', 'c2', m1=R1, m2=C2)

with monomers

Monomer('R1', ['x', 'c1'])
Monomer('R2', ['x', 'c1'])
Monomer('C1', ['y', 'c2'])
Monomer('C2', ['y', 'c2'])

The resulting ComponentSet is

ComponentSet([
Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=2) + R1(x=None) <> C2(y=1, c2=50) % C1(y=2) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr),
Parameter('bind_C2C1_R1_kf', 0.0001),
Parameter('bind_C2C1_R1_kr', 0.1),
Rule('bind_R1_C2', R1(x=None) + C2(c2=None) <> R1(x=1) % C2(c2=1), bind_R1_C2_kf, bind_R1_C2_kr),
Parameter('bind_R1_C2_kf', 0.0002),
Parameter('bind_R1_C2_kr', 0.2),
Rule('bind_R1R2_C2C1', R1(x=None, c1=1) % R2(x=1) + C2(y=2, c2=None) % C1(y=2) <> R1(x=50, c1=1) % R2(x=1) % C2(y=2, c2=50) % C1(y=2), bind_R1R2_C2C1_kf, bind_R1R2_C2C1_kr),
Parameter('bind_R1R2_C2C1_kf', 0.0003),
Parameter('bind_R1R2_C2C1_kr', 0.3),
])

The first rule in this set is erroneous as it has two dangling bonds:

Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=2) + R1(x=None) <> C2(y=1, c2=50) % C1(y=2) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr)

Looking into the code, the problem appears to be due to the use of the |= operator at line 825 of macros.py:

components |= bind_complex(s_row, row_site, s_col, col_site, klist, m1=m1, m2=m2)

Printing the progress to the screen, the output is:

----------Iteration 1----------
ComponentSet([
Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=1) + R1(x=None) <> C2(y=1, c2=50) % C1(y=1) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr),
Parameter('bind_C2C1_R1_kf', 0.0001),
Parameter('bind_C2C1_R1_kr', 0.1),
])
----------Iteration 2----------
ComponentSet([
Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=1) + R1(x=None) <> C2(y=1, c2=50) % C1(y=1) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr),
Parameter('bind_C2C1_R1_kf', 0.0001),
Parameter('bind_C2C1_R1_kr', 0.1),
Rule('bind_R1_C2', R1(x=None) + C2(c2=None) <> R1(x=1) % C2(c2=1), bind_R1_C2_kf, bind_R1_C2_kr),
Parameter('bind_R1_C2_kf', 0.0002),
Parameter('bind_R1_C2_kr', 0.2),
])
----------Iteration 3----------
ComponentSet([
Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=2) + R1(x=None) <> C2(y=1, c2=50) % C1(y=2) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr),
Parameter('bind_C2C1_R1_kf', 0.0001),
Parameter('bind_C2C1_R1_kr', 0.1),
Rule('bind_R1_C2', R1(x=None) + C2(c2=None) <> R1(x=1) % C2(c2=1), bind_R1_C2_kf, bind_R1_C2_kr),
Parameter('bind_R1_C2_kf', 0.0002),
Parameter('bind_R1_C2_kr', 0.2),
Rule('bind_R1R2_C2C1', R1(x=None, c1=1) % R2(x=1) + C2(y=2, c2=None) % C1(y=2) <> R1(x=50, c1=1) % R2(x=1) % C2(y=2, c2=50) % C1(y=2), bind_R1R2_C2C1_kf, bind_R1R2_C2C1_kr),
Parameter('bind_R1R2_C2C1_kf', 0.0003),
Parameter('bind_R1R2_C2C1_kr', 0.3),
])

Here we see that in the first two iterations the first rule is correct:

Rule('bind_C2C1_R1', C2(y=1, c2=None) % C1(y=1) + R1(x=None) <> C2(y=1, c2=50) % C1(y=1) % R1(x=50), bind_C2C1_R1_kf, bind_C2C1_R1_kr)

However, in the third iteration the C1(y=1) on both sides of the rule is replaced by C1(y=2). C1(y=2) appears in the third rule. Therefore, it appears that it is somehow overwriting the C1(y=1) in the first rule. I assume that this is due to the |= operator, although I'm not certain.

Ideally, BioNetGen would throw an error in cases like this. However, dangling bonds are currently allowed in BNGL. This behavior has been flagged for deprecation in a future release, however.

[alubbock]

This macro has been updated since Leonard wrote the issue. Is it fixed now? It no longer generates dangling bonds at least.

I've added dangling bond detection as a pull request here: pysb#328