A reaction rate that uses a Species observable stays frozen at its starting value

[wshlavacek]

I found a bug where NFsim does not update a reaction rate when that rate is
written as a function of a Species observable. When the same rate is written as a
function of a Molecules observable, it works fine. Here is a small model that demonstrates the behavior.

In the model below there is one molecule called A that starts at 100 copies and
slowly degrades. There are two products, Ps and Pm. Both are created at a rate equal
to the current number of A molecules. The only difference between them is that the
rate for Ps reads the A count through a Species observable, and the rate for Pm reads
the A count through a Molecules observable. Since both observables are counting the
exact same molecule, Ps and Pm should come out the same.

begin model
begin parameters
  k_deg 0.01
  k_prod 1.0
end parameters
begin molecule types
  A()
  Ps()
  Pm()
end molecule types
begin seed species
  A() 100
  Ps() 0
  Pm() 0
end seed species
begin observables
  Species   Sobs A()
  Molecules Mobs A()
  Molecules Ps_n Ps()
  Molecules Pm_n Pm()
end observables
begin functions
  fs() = k_prod*Sobs
  fm() = k_prod*Mobs
end functions
begin reaction rules
  Rdeg: A() -> 0  k_deg
  Rps:  0 -> Ps()  fs()
  Rpm:  0 -> Pm()  fm()
end reaction rules
end model

I made the XML for this with BioNetGen and then ran it with NFsim. The Species
observable requires the complex bookkeeping flag, so I ran with -cb and asked for a
single output line at the end:

NFsim -xml model.xml -sim 300 -oSteps 1 -o out.gdat -seed 1 -cb

Here is the result at the end of the run. The two observables Sobs and Mobs stay
equal to each other the whole time, but the two products are very different:

time   Sobs   Mobs    Ps_n     Pm_n
300      6      6      29885     9392

Ps came out at about 29885. That is what you get if the rate stays at its
starting value of 100 per second for the whole 300 seconds, because 100 times 300 is
30000. So the rate that reads the Species observable never noticed that A was
dropping from 100 down to about 6. It stayed frozen at the value it had at the start.
Pm came out at about 9392, which is the correct answer. I ran three different seeds and Ps was about 30000 every time while Pm was
around 9000 to 10000, so this is not just random noise.

One extra clue that might help you find it. If I ask NFsim to write more output lines
during the run, for example with -oSteps 4 instead of -oSteps 1, then Ps comes out
lower, around 13000 instead of 30000. It looks like the rate does get refreshed when
NFsim writes an output line, but not when a normal reaction changes the observable.
So writing output more often makes the Species version less wrong. A second variant where A instead grows from 0 up to 100 shows the same thing
from the other side.

I also ran the same model in RuleMonkey. There both Ps and Pm came out around 9200, so they agree
with each other and with the Molecules version from NFsim. That makes me
confident the model itself is fine and the problem is specific to how NFsim handles a
rate that depends on a Species observable.

Versions I used:
NFsim v1.14.3
BioNetGen 2.9.3
macOS