Setting really small values in ODE to zero

Specific Domains Modelling & Simulations
1

In modelling biological communities with ODEs, given that individuals are somewhat finite entities, I am trying to avoid infinitesimal densities which are mathematically correct but make no sense. More specifically, I want to set an extinction threshold at a certain small value (eg 1e-12), to ensure that extinct species are actually extinct and do not come back from the dead later. Currently doing this with a callback (see below). It works but I was wondering if there is a more elegant or efficient way of setting to zero values of the solution under a certain threshold?

My current approach (simplified, this is of course the drug concentration decay example with a large threshold)

function f(du,u,p,t)
    du[1] = -u[1]
    du[2] = -u[2]
end
threshold = 0.2
condition(u,t,integrator) = any( (u .< threshold) .& (u.>0))
function affect!(integrator)
  integrator.u[(integrator.u .< threshold) .& (integrator.u .> 0)] .= 0
end
cb = DiscreteCallback(condition,affect!)
u0 = [1, 0.5]
prob = ODEProblem(f,u0,(0.0,10.0))
sol = solve(prob,Tsit5(),callback=cb,saveat = 0.1)
plot(sol)
1 Like
2

I think this is a fine way to do this kind of modeling. Using Gillespie-type models are another way to have “true zeros” in the model.

3

Thank you for the reply :slight_smile: I tried discrete-individuals stochastic models (Gillespie algorithm, tau-leaping), but they need very large populations and very small mutation rates to approach the analytical solution. Also the added stochasticity was more an annoyance than anything, so I went back to ODEs.

4

makes sense