I saw that a mixed birth-death ODE-jump process I was using wasn’t behaving as expected.
Following some investigation, it appears that when I switch the rates from being ‘ConstantRateJump’ to ‘VariableRateJump’ the birth-death model stops properly simulating the theoretical expectation for a birth-death process.
This issue is also only noticeable when the population sizes are small, so I’m assuming that this means the issue is something to do with how the time between jumps is calculated when the gap between jumps is large.
Here is a simple toy birth-death jump process model I wrote to investigate the problem (the real problem is more complicated - so although in this version the variable rates aren’t necessary, they are in the full model. Also, in the full model, the rates do depend on part of the ODE solution, and hence the continuous problem has to be used).
function grow_fxn_sjm(n0::Float64, b::Float64, d::Float64, tmax::Float64)
u0 = [n0]
tspan = (0.0,tmax)
p = [b,d]
# ODE function
function ode_fxn(du, u, p, t)
b,d = p
du .= 0
nothing
end
function birth!(integrator)
integrator.u[1] += 1
nothing
end
function death!(integrator)
integrator.u[1] -= 1
nothing
end
b_rate(u, p, t) = (u[1] * p[1])
d_rate(u, p, t) = (u[1] * p[2])
b_jump = VariableRateJump(b_rate, birth!)
d_jump = VariableRateJump(d_rate, death!)
ode_prob = ODEProblem(ode_fxn, u0, tspan, p)
sjm_prob = JumpProblem(ode_prob, Direct(), b_jump, d_jump)
sol = solve(sjm_prob, Tsit5())
return sol
end
If I run this model many times with a starting population (n0) of 1.0 , I can see that the model is under-estimating the expected population size at time t, given some birth and death rates b & d. For the simple birth-death process, this should simply be n(t) = n0 * e^((b - d)*t) - if I change the b_jump and d_jump to ConstantRateJumps then the model does on average reach the expected population size, (nt).
Does anybody know if I’m not using the VariableRateJumps correctly? I’ve read the documentation on using upper and lower bounds and a dependency graph, but I think because I need the problem to be continuous I can’t use these features.
Any help much appreciated!
