gillespie doesn't update time-dependent assignmentRules in reaction rates

[wshlavacek]

When I run BIOMD0000001040 with gillespie SSA, the mean trajectory doesn't match the ODE solution
even though it probably should. It looks like an assignment rule that depends on time is not being re-evaluated
as the simulation advances, so the propensity it feeds is stale by one output interval.

The model has an assignment rule Mpl = Aexp(ct) - Bexp(dt), which is an explicit function of time,
and Mpl appears in the rate of the synthesis reaction that produces Mrbc. With CVODE this behaves
correctly. With gillespie, Mrbc stays frozen at zero for the first output interval and then trails the
ODE for the rest of the run. I am on roadrunner 2.9.2.

import numpy as np, roadrunner
rr = roadrunner.RoadRunner("BIOMD0000001040.xml")

ODE

rr.timeCourseSelections = ["time", "[Mrbc]"]
print("ODE:", np.array(rr.simulate(0, 10, 6))[:, 1])

gillespie, mean over 40 replicates

means = []
for s in range(40):
rr.reset()
rr.setIntegrator("gillespie")
rr.integrator.seed = s + 1
rr.timeCourseSelections = ["time", "[Mrbc]"]
means.append(np.array(rr.simulate(0, 10, 6))[:, 1])
print("SSA:", np.mean(means, 0))

This is what I get:

ODE: [0, 36.53, 83.08, 113.69, 129.62, 135.46]
SSA: [0, 0, 49.67, 83.12, 101.67, 109.90]

For this model Mrbc is effectively an immigration–death process with a time-varying immigration rate,
so the SSA mean has to converge to the ODE solution exactly, and CVODE confirms that the ODE values
above are the correct answer. The gillespie mean is shifted down by roughly one output interval, which
is consistent with the propensity being computed from Mpl evaluated at the previous sample time rather
than the current time.

[luciansmith]

Is this not a duplicate of #1317?

[wshlavacek]

Yes — I'd treat #1318 as a minimal reproducer of the same underlying bug as #1317, so consolidating them makes sense.

The common defect is that gillespie doesn't re-evaluate time-dependent assignmentRules when it computes reaction propensities, so a rate that depends on such a rule lags by one output interval. #1318 isolates exactly that with a single explicit function of time (Mpl = A·exp(c·t) − B·exp(d·t) feeding the Mrbc synthesis rate): CVODE grows Mrbc from the start, gillespie stays frozen at 0 for the first interval and then trails for the rest of the run. For that affine immigration–death system the SSA mean must equal the ODE exactly, so the gap is unambiguous.

#1317 is the original, harder-to-read report — a BNG-generated 10-species model where the same stale-rule mechanism is tangled with other features (I couldn't isolate which), surfacing as three different biases on three by-construction-identical toggles rather than a clean one-interval lag. That non-uniform pattern is actually what you'd expect from a stale time-dependent rate sampled at different stochastic times, not a uniform error like a missing volume factor.

Suggestion: keep #1318 as the canonical reproducer to develop and test the fix against — it has an exact analytical oracle — and confirm the fix also clears #1317.