Kappa-Dev / KappaTools

Tool suite for kappa models. Documentation and binaries can be found in the release section. Try it online at
http://kappalanguage.org/
GNU Lesser General Public License v3.0
112 stars 41 forks source link

Unexpected positive self-influence in DIN with explicit unary rates #592

Open slegare2 opened 5 years ago

slegare2 commented 5 years ago

Hi,

Binding rules show in unexpected positive self-influence in DIN when combining explicit unary rates and branching polymers.

The following example contains a single agent type A which can bind to itself on 3 different sites. It can hence form large branching polymer. Each binding rule has an explicit unary rate {0}.

//
//     A
//     |
//     A
//    / \
//   A   A
//

%agent: A(x y z)

%init: 100 A()

'Ax binds Ax'       A(x[./1]), A(x[./1]) @ 0.01 {0}
'Ay binds Ay'       A(y[./1]), A(y[./1]) @ 0.01 {0}
'Az binds Az'       A(z[./1]), A(z[./1]) @ 0.01 {0}

'Ax unbinds Ax'     A(x[1/.]), A(x[1/.]) @ 0.01
'Ay unbinds Ay'     A(y[1/.]), A(y[1/.]) @ 0.01
'Az unbinds Az'     A(z[1/.]), A(z[1/.]) @ 0.01

%obs: 'AxAx'     |A(x[1]), A(x[1])|
%obs: 'AyAy'     |A(y[1]), A(y[1])|
%obs: 'AzAz'     |A(z[1]), A(z[1])|

%mod: alarm 0.0001 do $DIN "flux_triple.json" [true];
%mod: alarm 100 do $DIN "flux_triple.json" [false];

When I run this model in the Kappapp, I get the following DIN. We see the positive self-influence on all 3 binding rules. Note that the static influence map presents no such positive self-influences.

flux_pos

However, if I lower the binary rate of each binding rule to 0.0001, the expected negative self-influence is recovered (at least qualitatively).

flux_normal

The previous positive self-influences seem like a bug to me. Although I do not know very well how KaSim works under the hood, I can hypothesize that it tries to bind two As, and then breaks the bond when it realizes they are part of the same complex. It would then somehow count this as a bond break, hence the creation of a new possible bond to form. Something like that?

Thanks, Sebastien

pirbo commented 4 years ago

Here is what happen, The way the simulator works with respect to molecular ambiguity is by over-approximating most of the time the number of binary instances by the total number of instances while computing the activity of the mixture but by discarding dynamically selected an instance that happen to be unary (doing what is called a null event). BUT, if too many null events appears in a row then the exact computation of the number of binary instances is computed (that decreases the computed activity of the binary rule which now the exact one). Then, the first time that rule is applied afterward, its activity is recomputed (as it should) but we are back into computing the over-approximated activity (which is bigger then the former exact activity you had) and the DIN mechanism ignores this subtlety and store delta_activity = new_over-approximated_activity - old_exact_activity. That's how you end up with a positive number in the delta_activity a rule has on itself here...

slegare2 commented 4 years ago

Thank you for looking into this! Could this be resolved by still computing the over-approximated activity when the exact computation is done and then doing delta_activity = new_over-approximated_activity - old_over-approximated_activity ?

pirbo commented 4 years ago

I have a plan B which is still incorrect but maybe less incorrect:

If I register in the DIN the delta_activity "exact_activity - over-approximated_activity" when I compute the exact activity it fixes your case.

The case that becomes incorrect is the one where the rule is never applied after having been refined because there will be a malicious self-influence that in fact correspond to the refinements...

pirbo commented 4 years ago

OK ... Sight ... of course ... Sorry ....

If in "absolute" mode exact_activity - old_approx_activity + new_approx_activity - exact_activity is what you expect, in "relative" mode (the default) (exact_activity - old_approx_activity)/old_approx_activity + (new_approx_activity - exact_activity)/exact_activity is ... pure non sense !