A bug?

[Vilin97]

I test the network A+B <--> C with forward rate 0.1 and backward rate 1 to time 10 seconds, starting at 500 A and B and 0 C. I consistently (average of 10^5 simulations) get that in the end there are ~0.12-0.29 less A and B molecules than there should be via an analytic equilibrium calculation, which is 65.887. If it were the case that the number of A and B molecules is larger than in steady state, I would conclude that the network just does not fully reach the equilibrium. But the number of A and B is lower than in equilibrium, meaning the simulation "overshoots" the equilibrium. 0.22 error is within 0.4% relative tolerance, which seems low but the consistency of this result is worrisome. I see this behavior with all SSAs I tested: RSSA, RSSACR, NRM, SortingDirect, RDirect. Is it possible that there is a statistical bug somewhere? It would have to be in one of the common parts of SSAs since the behavior is consistent across different SSAs.

[isaacsas]

Can you post your code? How many simulations did you run? (Missed you said 10^5 simulations.)

[isaacsas]

What happens if you run to a larger final time?

[isaacsas]

What is your analytic mean formula? The formula from McQuarrie (1967) agrees with the formula in "Deterministic Versus Stochastic Modelling in Biochemistry" after page 74 (though they could have taken it from his paper). One caveat; I think there is a typo in McQuarrie and K should be k2/k1 not k1/k2. This is also consistent with the book's table. I get the mean should be 65.6708109088797.

[isaacsas]

We should add this as a test anyways...

[isaacsas]

using Catalyst, DiffEqJump, HypergeometricFunctions

Nsims = 1e5

rn = @reaction_network begin
    .1, A + B --> C
    1.0, C --> A + B
end
u0 = [500,500,0]
tspan = (0.0,10.0)
dprob = DiscreteProblem(rn,u0,tspan)
jprob = JumpProblem(rn,dprob,Direct(),save_positions=(false,false))

function getmean(jprob,Nsims)
    Amean = 0
    for i = 1:Nsims
        (mod(i,div(Nsims,100)) == 0) && println("i = $i")
        sol = solve(jprob,SSAStepper())
        Amean += sol[1,end]
    end
    Amean /= Nsims
    Amean
end

Amean = getmean(jprob, Nsims)

K = 1/.1
function analyticmean(u, K)
    α = u[1]; β = u[2]; γ = u[3]
    @assert β ≥ α "A(0) must not exceed B(0)"
    K * (α+γ)/(β-α+1) * pFq([-α-γ+1], [β-α+2], -K) / pFq([-α-γ], [β-α+1], -K)
end
analytic = analyticmean(u0, K)

println("Amean = $Amean, analytic = $analytic")

[Vilin97]

65.6708109088797 is very close to the mean I get from simulations! Closing the issue