Zombies in biological ODE : why is my solver not sticking to zero?

Specific Domains Numerics
, , ,
1

Hi numerical ODE solvers enthusiasts, we’ve got a problem here that we think necessitates your kind input :slight_smile: Here are two related questions:

The following ODE system represents biological species biomasses B, and it is theoretically guaranteed that species with null biomass B = 0 have a constant, null trajectory.
However, and depending on how we parametrize DifferentialEquations’s solver, we get “zombie” species rising out from the dead to nonzero biomasses even if their initial state is zero. We are not sure why.

using DifferentialEquations

"""
ODE system defining the evolution of the biomasses of my 3-species system.
Species 2 and 3 feed on species 1.
"""
function dBdt!(dB, B, p, t)
    println("t: $t, B: $B") # <- Sometimes called with B[2] > 0 : why?
    # @assert B[2] == 0 # <- This would fail, why?
    r, K, B0, x, y2, y3 = p # unpack parameters

    # Compute feeding rates
    F21 = (x * y2 * B[2] * B[1]^2) / (B0^2 + B[1]^2)
    F31 = (x * y3 * B[3] * B[1]^2) / (B0^2 + B[1]^2)

    # Compute growth rates
    dB[1] = r * B[1] * (1 - B[1] / K) - F21 - F31
    dB[2] = F21 - x * B[2]
    dB[3] = F31 - x * B[3]

    # We have (B[2] = 0) ⇒ (dB[2] = 0),
    # but even (B[2] = 0) is not always true, why?

end

# Second species biomass is set to zero, expected to stay at zero.
B0 = [0.1, 0.0, 0.2]
p = (1, 1, 0.5, 0.314, 7.992, 8)
problem = ODEProblem(dBdt!, B0, (0, 100_000), p)
no_zombies = solve(problem, reltol=1e-6) # All B[2] values are zero.
zombies! = solve(problem, reltol=1e-3) # Some B[2] values are becoming positive.

Sneaking into the above example, we realized that dBdt! is sometimes called with non-zero values of B[2]. As a consequence, guaranteeing that B[2] = 0 ⇒ dB[2] = 0 is not enough to prevent the apparition of zombies.

Here is my first question : Why would solve ever call dBdt! with non-zero values of B[2] in the above example?

Now, here is more general context: dBdt! is actually generated from user input (foodweb, species capabilities, etc.). As a consequence, even though zombies did not appear in the particular above example when reltol = 1e-6, there is no way we can calculate generic, adequate values for reltol upfront.

In fact, the real need we have is exactly the one described in this other thread : Given a positive extinction threshold xt > 0, trigger a callback as soon as B[i] < xt is detected and force it to B[i] = 0 instead… and forever.

In this other thread, @ChrisRackauckas suggested that using callbacks from DiffEqCallbacks to this end was a fine approach. However, the example above demonstrates that hard-setting B[i] = 0 is not a guarantee that the trajectory of B[i] will stick to 0 forever.

So here is my second question : Are DiffEqCallbacks really the right tool for the job?
If yes, then how to use them properly and avoid zombie species?
If not, should we rely on other approaches, like stopping simulation whenever B[i] < 0 is detected then re-generating dBdt! with fewer variables?

1 Like
2

How about altering the equations so for small values of B it’s attracted to 0? This is more of a realistic model anyway. You probably need some nonlinear term like -B(1-inverse_logit((B-B0)/s)) whose effect goes away for larger B but at small B has a downward drag

2 Likes
3

Another thing you can try is using a callback set B to exactly zero when it falls below some threshold. You can modify states in a callback through the integrator interface

4

Nice idea. Just beware that the narrower the gap of that canyon this introduces around B=0, the stiffer the equations of motion will be.

I think the best answer to this question is either not be bothered by slow, small creep away from B=0, or to do a discrete simulation, where populations take on only integer values, instead of an ODE.

It seems like if your ODE model is realistic, the eigenvalues of the linearized growth rates of each species near zero would be comparable in magnitude, meaning it’d take a relatively long time, compared to other species’ life cycles, for floating-point or time-stepping errors to grow many orders of magnitude to appreciable size.

1 Like
5

I’m a huge fan of agent/individual based simulations especially in biological applications. This is a good idea for population dynamics.

1 Like
6

BTW For those who are interested in zombies in ODEs here’s an oldie from my blog back in the day… the LaTeX plugin I was using stopped being maintained… unfortunately.

http://models.street-artists.org/2010/03/01/improved-zombie-dynamics/

7

And, presciently… here’s an agent based zombie model I built with my kids for Christmas 2019… Immediately before another less well known pandemic…

http://models.street-artists.org/2019/12/25/merry-zombie-apocalypse/

9

This is what @iago-lito discuss at the end of the post. Even with such a callback (that we already tried) there is no guarantee that a species set extinct with a callback B[i] = 0.0 will stick to 0 as illustrated in the example above

1 Like
10

Thanks for the suggestion but we are currently developing a package for a well-known and commonly used model in ecology: the bio-energetic model (see for instance Yodzis 1992 or Williams, Brose & Martinez 2007). Thus, it isn’t an option for us to switch for agent based models.

Concerning your suggestion of adding a term to the equation, first as said above we are writing a package implementing a widely used model so we don’t have room to modify it (or very little). Anyway, I think this misses the core question which is why dBdt! is fed with non-zero value whereas the biomass is exactly 0? We would like to understand what is happening there (in solve) and how to deal with it properly. That is why saying don’t care about ‘small creep away of B=0’ does not seem very satisfying. In the minimal example it has very little effect, but we start to run complex simulations with 100+ species, it is not guarantee that small creeps won’t have any significant effect on the other species dynamics.
Also, as a side-note, the example given in @iago-lito post is an extreme case where the species extinction happens very slowly and we use a very long time span which goes from t=0 to t=100_000, consequently the solver takes very large time-steps. However, most of the time species go extinct way faster. But still we have to be ready for such cases.

1 Like
11

Hi @John_Gibson and thank you for feedback. Are you meaning that the small positive numbers that we see creeping up in trajectories supposed to be null would never actually become big and always stay close to zero? If yes, this would be very good news for us, but what guarantee do we have?

I totally agree with you that discrete-population models are also very-well suited to represent biological populations. However, they cannot beat ODEs regarding performance when many generations are needed, or mathematical analysis opportunities for small, tractable foodwebs. Also, our current project is to implement the particular ODE models cited by @ismael-lajaaiti above, which are widely used in our community. So the problem is not about picking the adequate representation (ODE vs. MAS), but about correctly implementing this particular ODE representation.

Hi @dlakelan and thank you for feedback as well. I like your idea very much, although it would complicate the functional form of our ODE systems in a way that possibly makes it impossible to get analytical results for small foodweb modules.

This being said, I can think of another “hybrid” representation, very close to yours, where the population biomass is continuous, but the simulation domain is only D_1 = \{0\}\ \cup\ [x, +\infty[ with x>0 instead of D_2 = \mathbb{R}^+. Although our strict ODEs do describe evolution of trajectories within D_2, we can restrict them to D_1 by deciding that, if any argmin_t(B[i](t) < x) exists, then further trajectories are better described by a new ODE system with variable B[i] removed and B[i](t) = 0\ \forall t \in \mathbb{R}^+.

The above model preserves the structure of the existing, consensual ODE models while still making B “attracted to 0” like you said, and yet without making the ODEs stiffer like @John_Gibson was worried about.
Furthermore, its implementation is simple as it just boils down to stopping the solver whenever B[i] < x is detected, then restarting it with updated equations.

Our (unanswered) concern is that we initially thought that a more simple implementation of the above was possible, using DiffeqCallbacks to force B[i] = 0 without having to restart the solver. Is it? Or would it be wrong? And why does solve call our dBdt! function with non-zero values of B[2] in the above example?

2 Likes
12

Ok, so this is useful information. you’re trying to implement someone elses equations, for all sorts of valid reasons, and you want to use the “Callback method” to just extinctify a species if it drops below some threshold value… reasonable.

I don’t know why the solver calls your dBdt! function with nonzero B[2] but it seems like a callback should allow you to extinctify the B[2] population whenever it hits some threshold, let’s say for example 1e-8 on a dimensionless scale where a typical population is 1. However, i’m not at all convinced that this won’t make the system “stiff” in the general sense that to determine whether the callback should be called, the system will have to track the B population to 1 part in 1e8 and it may need to take short timesteps or do other extra calculations to ensure that your population doesn’t just “jump” over that threshold and continue to grow.

I don’t know enough about DifferentialEquations.jl to be honest. It’s one of those parts of Julia I am itching to work with, but haven’t had a good reason to do so yet. The one project I wanted to work on turned out to be a bit too ambitious (gravitation with thousands of stars and retarded gravity, so a DAE with millions of couplings).

1 Like
13

it may not stick to zero, but it will constantly be forced to zero. If you trigger the callback at a very small value then you can bound the maximum error in B[i] but if your equations are sensitively dependent on any of the B[i] having a nonzero value even for a small time may produce an error in the other populations that grows exponentially. An example might be a population of organisms that produce a potent long-lasting toxin. Other populations may plummet if B[n] goes above zero even for a short time.

1 Like
14

If you try using an explicit solver like Tsit5 this issue goes away:

zombies! = solve(problem, Tsit5(), reltol=1e-3)

Maybe the linear or nonlinear solve embedded in an implicit method is what is introducing the small perturbations away from zero?

1 Like
15

Switching the linear solver to a Krylov method seems to avoid the issue for implicit methods (in limited testing):

using LinearSolver
zombies! = solve(problem, Rosenbrock23(linsolve=IterativeSolversJL_GMRES()), reltol=1e-3)
16

So I’ve got an interesting related problem. I coded up my zombies model:


using DifferentialEquations, StatsPlots

function zombies!(du, u, p, t)
    (H,Z,D,G,hwin) = u
    (b,d,r,e,rot,l,ldec,teach) = p
    du[1] = b*H - d*H - (1-hwin)*H*Z
    du[2] = r*D - hwin*H*Z + (1-hwin)H*Z -e*H*sqrt(Z)
    du[3] = d*H - rot*D - r*D
    du[4] = hwin*H*Z + rot*D + e*H*sqrt(Z)
    du[5] = l * hwin*sqrt(1-hwin)*H*Z - ldec*hwin+teach*(1-hwin)
end

If I solve it with some “base” values, it gives reasonable solution with sufficiently high tolerance:


zinit = 1e-5
hwininit = .1

# We simulate time in days because its makes estimating some of the rates easier
pvals = [.065/365,.02/365,.001/365,0.00,.15,5.0,.1,0.0]
prob = ODEProblem(zombies!,[1.0-zinit,zinit,0.015,0.0,hwininit], 60.0, pvals)

try
    sol = solve(prob,Tsit5(); abstol = 1e-11, reltol=1e-8)
catch e
    println(e)
end

plot(sol; title="Humans u[1] vs Zombies u[2]")

image
image849×549 29.7 KB

If I try to rescue the population by adding in my zombie killer elite: e > 0

pvals2 = copy(pvals)
pvals2[4] = .05
prob = ODEProblem(zombies!,[1.0-zinit,zinit,0.015,0.0,hwininit], 60.0, pvals2)


try 
   sol = solve(prob,Tsit5(); abstol=1e-11,reltol=1e-8)
catch e
    println(e)
end
plot(sol)

It borks when Z drops below 0.0 with a sqrt of negative number error.

If I try a callback that sets Z to 0.0 when Z drops below 1e-6


zombieextinct = ContinuousCallback((u,t,int) -> u[2] - 1e-6,
    nothing,
    int -> int.u[2] = 0.0)

try
    sol = solve(prob,Tsit5(); 
        callback = zombieextinct, abstol=1e-11, reltol=1e-8)
catch e
    println(e)
end

plot(sol)

It still tries to take the sqrt(Z) when Z is negative.

Note that D >=0 and H >= 0 implies that at Z=0 dZ/dt >=0 so it should never go negative.

Any suggestions?

17

I may have misunderstood your question, but if it is about forcing a state variable or species to 0 at all times or most of the time, I can think of three additional ways of making species at 0 if possible.

  1. Append the initial conditions as additional parameters to p. Then, you just include an if statement in the model that states that if p[inital conditions] == 0., then dB[your choice]/dt = 0. You can also propagate zeros to other equations by setting B[your choice] = 0.

  2. Instead of

dB[2] = (x * y2 * B[2] * B[1]^2) / (B0^2 + B[1]^2) - x * B[2]

, you can do

dB[2] = B[2]*((x * y2 * B[1]^2) / (B0^2 + B[1]^2) - x)

. Notice how I multiplied everything by B[2]. You can combine 1. with 2. to ensure dB[2] is always returning a zero value. This is the easiest way.

  1. Adjust absolute tolerance to a very small value. See Absolute and relative tolerance definitions - MATLAB Answers - MATLAB Central for more details.

I don’t think B[2] will ever be exactly 0 even if it shows as 0 (comparing floating point to integer doesn’t turn out the way you expect it) because of how floating point works. Anyways, combining the methods (or using them individually) above will significantly reduce the amount of nonzeros you see at the end of the simulation with your println statement. You don’t have to use all three. Hope this helps!

If you want to know why the solve sometimes spits up nonzero values, I think you can play around with the integrator interface instead of solve.

Note: Depending on how you define zeros, you may get type instabilities which affects the performance of the code.

18

I didn’t try this, but you can use an if statement before du[2] to set Z to zero.

19

yeah, that’s what I’m trying right now, as well as a callback for if Z < 1e-6 setting it immediately to 0.0

20

Hi @-all and thank you again for your feedback :slight_smile: Let’s attempt to sort all this out ^ ^"

@isaacsas Indeed, and also @kcqpo if we tweak reltol, abstol and/or other parameters of solve (callbacks aside). Unfortunately, we cannot just be satisfied with these approaches because the problem is not about avoiding zombies this one particular ODE system. We are willing to avoid them in any biological system possibly generated from user input. (Note that: no matter the user input, they all have this property that B[i] = 0 \implies dB[i] = 0, and it is our responsibility to enforce it.)

Since we do not currently understand why some parametrizations of solve() yield zombies and some do not, we have no guarantee yet that forcing Tsit5 or Rosenbrock23 or even abstol=1e-22 will always protect users against zombies. In addition, although providing good default parametrization for solve is a must, we would prefer not to force it so our users remain ultimately able to parametrize the solver the way they prefer.

Put it another way: There is no parametrization of solve() (abstol, reltol, alg, etc. callbacks aside) where zombies make sense. So we think we should not be avoiding zombies by relying on one particular parametrization rather than another.

@kcqpo This could work indeed, at the cost of one extra branch within dBdt!. I’m not exactly worried about performance yet, because if is_extinct(p, i) is easily predictable. However, and as it turns out, we (sometimes) have no data in p and the whole dBdt! function itself is generated from p using Julia metaprogramming features (effectively improving performances by x10 to x500 depending on the number of species). As a consequence, if we were to change the value of p during the simulation to extinguish species i, then we’re better off re-generating dBdt! with the whole variable dB[i] removed :wink:

@kcqpo This is not a guarantee that dB[2] be zero because, as it turns out, dBdt! gets sometimes called by solve() with non-zero values of B[2]. This is the very point of this topic: we do not understand why. Why would it “probe” our derivative function with non-zero values for a variable whose value and derivative have always been zero since the initial state? :frowning: Any clue appreciated ^ ^"

@dlakelan He-hee ^ ^ In the system you describe, “zombies” are actually desired but they mean a different thing. This being said, I think I can see why the problem described here could avoid you from having undesired populations crawling up from zero (or “meta-zombies”) :wink:

@dlakelan Very true. This is why the problem described here is of high relevance to us. B[2] = 0 and B[2] = 1e-6 have very different biological meaning and potentially lead to very different outcomes. We cannot tolerate zombies.

@dlakelan Sigh, no worries <3 :') Maybe someone else here does? Or else we maybe need to get deep into the solving algorithms internals/source code to understand what’s at play here? Could it maybe even be a… bug? :frowning:

@dlakelan This is very interesting and I didn’t know we have to think about stiffness this way. FWIU now, using a callback to force this “absorbing” behaviour of the value 0 is actually a creepy way of changing what dBdt! actually means, and making the ODE system actually simulated different from the system described by dBdt! alone.
Again, I think this is a good reason to embrace that our biological model is better represented by a sequence of (B0, dBdt!) pairs rather than a single one. When a species goes extinct, we maybe should just switch to another ODE system.

21

Hey @iago-lito, I took a deeper dive into @isaacsas solution. I am not a full expert on this by any means and correct me if I’m wrong guys.

Long Explanation: I think the reason why “zombies” appear when the initial condition is set to zero is because of the :auto option (see Common Solver Options (Solve Keyword Arguments) · DifferentialEquations.jl). If :auto is used, the backend chooses either
AutoTsit5(TRBDF2(autodiff = false)) or AutoTsit5(Rosenbrock23(autodiff = false)) as the solver. You can look at DifferentialEquations.jl/ode_default_alg.jl at master · SciML/DifferentialEquations.jl · GitHub for more details about the default choice. Basically, the algorithm switches between a stiff implicit solver and a nonstiff explicit solver. When the solver is an implicit method like Rosenbrock23, the algorithm is basically trying to find a root to a function. The root is the next state. For the optimization algorithm, it can be something like Newton’s method. The optimization algorithm will make many calls/guesses to your provided ODE function, and it does not guarantee 0.0 as the solution for the root. If some stiffness is detected, the algorithm switches from explicit to implicit. As a result, nonzero values appear when tolerances are low. See more at Solving Stiff Ordinary Differential Equations - MIT Parallel Computing and Scientific Machine Learning (SciML).

To verify my hypothesis (issue comes from the implicit solver), I compared an explicit method like Tsit5() to an implicit method like Rosenbrock23() for low tolerances (i.e. 1e-2, 1e-3, 1e-4). I also changed the assertion statement from @assert B[2] == 0 to @assert B[2] == 0.0 because we are working with floating point numbers. The result confirms my hypothesis that the issue stems from using implicit solvers. Using explicit solvers such as Tsit5() and RK4() do not assert an error, but implicit solvers Rosenbrock23(), Rodas4(), and KenCarp47() do. The implicit solver is correctly implemented as it should, but due to the optimization used, an exact integer solution like zero is most likely not going to happen if the tolerance is too low.

Short Answer: Based on the formulas for explicit solvers (where every function call is based on the previous state), I am confident that they should always result in zeroes for B[2] in your system. Use explicit solvers to guarantee a zero solution whenever the initial value is zero for your ODE problems (provided that every term is multiplied by zero). Implicit solvers uses optimization like Newton’s method to guess your next state. Optimization never guarantees an exactly zero solution. If implicit solvers must be used to tackle stiffness, then the problem won’t be avoided easily. You can try to solve stiff problems with explicit solvers, but it will take very long or sometimes fail if the problem is very ill-conditioned.

Hope this answers your first question about zombies when the initial condition B_{2, t=0}=0!

4 Likes