alubbock
commented
6 years ago 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/pysb#328
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.