Gillespie gives the wrong mean when a reaction rate depends explicitly on time.

[wshlavacek]

Fixes #1318. The direct method samples the waiting time from a propensity it treats as constant until the next reaction, which is only correct when the rate laws are time-homogeneous. When a rate depends explicitly on time (directly, or through a time-dependent assignment or rate rule), the SSA mean drifts off the ODE: a rate that is zero at t=0 froze the first reporting interval, and coarse output intervals trailed the true mean by about one interval. Now we check once, on the first step, whether any rate varies with time, and if so we find the next-reaction time by integrating the propensity over the wait instead of freezing it. Time-homogeneous models take the original code path unchanged.

Symptom

In BIOMD0000001040, Mrbc is produced at a rate proportional to Mpl(t) = A*exp(c*t) - B*exp(d*t), an explicit function of time. CVODE reproduces the ODE; gillespie stayed at 0 through the first output interval and trailed afterward.

simulate(0, 10, 6), mean of 60 fixed-step replicates:

ODE :  0   36.5   83.1   113.7   129.6   135.5
before: 0    0     49.7    83.1   101.7   109.9
after:  0   36.5   83.1   113.7   129.6   135.5

How it works

• On the first integrate call, compare each reaction rate at the current state for a few different times. If none change, the model is time-homogeneous and the original direct method runs (this detection adds no random draws, so those results are bit-identical to before).
• If a rate does vary with time, the waiting time tau satisfies integral over [t, t+tau] of s(u) du = -log(r1), where s is the total propensity re-evaluated as time advances (species are frozen between reactions). We accumulate that integral with an adaptive trapezoidal sub-step, then select the reaction using the propensities at the firing time. This reduces exactly to tau = -log(r1)/s when s is constant.

This is the direct-method form of the integrated-propensity representation in Anderson (2007): reaction times are the firings of unit-rate Poisson processes whose internal clock is the integral of the propensity. That representation is exact. This implementation realizes it by accumulating the integral with the trapezoidal sub-step, so the realized sampler is exact in the limit of small steps and converges to it as the tolerance is tightened; it is not bit-exact sampling at a finite tolerance.

Tests

• New regression test TimeDependentPropensityMatchesODE: an immigration-death process with a time-dependent immigration rate, whose ensemble mean has a closed-form ODE solution; the fixed-step SSA mean is checked against it.
• The time-homogeneous path is byte-for-byte identical to develop (verified on a mass-action model over many seeds).
• The stochastic test suite is unaffected: every model in it has time-homogeneous rate laws (the CSymbolTime cases use time only in event triggers), so they all take the unchanged path. Spot-checked models 00028 and 00029 produce output identical to develop.

Notes

• The sub-step is bounded by a relative-change tolerance on the propensity, a fixed constant (0.05) in the code. Accuracy improves as it shrinks, and it has no effect on time-homogeneous models.
• The time-dependence check evaluates each rate at the same state for a few different times and looks for a change. That is a value-based heuristic, and can in principle miss a rate that is flat at the sampled times but varies elsewhere. Inspecting the rate-law expression trees for the time symbol (following the assignment-rule graph) would be more robust than sampling values; the code notes this.
• Events are evaluated at reaction times, as before, so an event timed inside a long quiet interval fires at the next reaction, the same granularity as the existing solver.
• This is independent of gillespie simulator: wrong stationary distribution on a 10-species SBML with piecewise/assignmentRule kinetic laws and stochastic two-state toggles #1317 (the other Gillespie fix) and can be merged in either order.

Reference

David F. Anderson, "A modified next reaction method for simulating chemical systems with time dependent propensities and delays," J. Chem. Phys. 127, 214107 (2007). arXiv:0708.0370

[luciansmith]

This is interesting, because for everything else, the gillespie simulator just claims that it cannot handle elements changing in time--models with rate rules cannot be simulated. As you point out, we don't cover the term 'time', so if 'time' appears in a reaction rate or in an assignment rule, the result will be incorrect.

My instinct would normally be to just look for 'time' and treat it as an error in a gillespie simulation, if it showed up in assignment rules or reaction rates. But we already have this fix for 'time' in reaction rates--@hsauro, do you think it's worth it to keep this, and just error out for 'time' in assignment rules?

(Alternatively, I suppose we could try to handle 'time' in assignment rules, too, but that seems difficult.)

[luciansmith]

Wait, some of the comments in the code seem to indicate the time-dependent assignment and rate rules would be handled by this change. There's a lot of math, so I don't pretend to understand it, but is this true?

[wshlavacek]

The code handles explicit time dependence in a reaction rate when the time dependence enters directly or through an assignment rule. Rate rules are not handled. Sorry for the confusion. The misstatement has been fixed in the code comments.

[wshlavacek]

I agree that an error would be better than silence. That said, with this change, the SSA now yields results consistent with the ODE solver.

[wshlavacek]

If you want rate rule support, it's a natural extension of the Anderson method implemented here. Advance the rate rule variable through the waiting interval via getStateVectorRate used by CVODE. It's about 100 lines plus the trouble of testing.