Help with misbehaving ODE simulation

Specific Domains Modelling & Simulations
,
1

Hi all,

I am having some issues with a set of ODEs I am simulating with DifferentialEquations.jl and am wondering if anyone can help. Basically, the code is simulating populations of bacteria and their viruses (phage) that are interacting, with interactions depending on certain characteristics (d_elems, cd_elems). Half of the bacterial strains only have these characteristics when total bacterial abundance is above a threshold, the others always have them. Attached is the code I am using, but the solution seems to be unstable for certain values of the threshold (such as 1e10, the value in the code currently), even though (I think) the equations in principle should not be ill-behaved. Any suggestions on what might be going on (is there a different solver I should be using?) would be incredibly helpful. I can provide more context for the equations if needed. Thanks!

using DifferentialEquations, Combinatorics

function generate_subsets(n, size)
    return collect(combinations(1:n, size))
end

function infection_rate(k, subset1, subset2)
    return (subset1 ⊆ subset2) * k
end

function model(du, u, p, t)
    δ, k, µ, λ, ν, b_max, α_max, n_max_b, n_max_p, n_bacteria, n_phage, d_elems, cd_elems,
    threshold = p
    B = u[1:n_bacteria]
    P = u[n_bacteria+1:end]

    active = sum(B) >= threshold

    g(nb) = α_max * cos(π * nb / (2 * n_max_b))
    h(np) = b_max * cos(π * np / (2 * n_max_p))

    for (i, subset_b) in enumerate(d_elems)
        αi = g(length(subset_b))
        infection_sum = sum(infection_rate(k, subset_b, subset_p) * P[j] for (j, subset_p) in enumerate(cd_elems))
        du[i] = B[i] * (αi - µ - infection_sum) + λ
    end

    for (i, subset_b) in enumerate(d_elems)
        αi = active ? g(length(subset_b)) : α_max
        subset_b = active ? subset_b : []
        infection_sum = sum(infection_rate(k, subset_b, subset_p) * P[j] for (j, subset_p) in enumerate(cd_elems))
        du[i + length(d_elems)] = B[i + length(d_elems)] * (αi - µ - infection_sum) + λ
    end

    for (j, subset_p) in enumerate(cd_elems)
        bj = h(length(subset_p)) - 1
        mutant_infection = active ? sum(infection_rate(k, subset_b, subset_p) * B[i] for (i, subset_b) in enumerate(d_elems)) : sum(infection_rate(k, [], subset_p) * B[i] for (i, subset_b) in enumerate(d_elems)) 
        infection_sum = sum(infection_rate(k, subset_b, subset_p) * B[i] for (i, subset_b) in enumerate(d_elems)) + mutant_infection 
        du[n_bacteria + j] = P[j] * (-δ + bj * infection_sum) + ν
    end
end

function main()
    δ, k, µ, λ, ν, b_max, α_max = 1.0, 1.0, 1e-2, 1e-15, 1e-15, 15, 10.0
    ntot, n_max_b, n_max_p = 6, 3.27, 4.5

    d_elems = generate_subsets(ntot, 2)
    cd_elems = generate_subsets(ntot, 3)

    n_bacteria = length(d_elems) * 2 
    n_phage = length(cd_elems)
    threshold = 1e10 

    p = (δ, k, µ, λ, ν, b_max, α_max, n_max_b, n_max_p, n_bacteria, n_phage, d_elems, cd_elems, threshold)
    u0 = ones(n_bacteria + n_phage)
    tspan = (0.0, 100.0)

    prob = ODEProblem(model, u0, tspan, p)
    sol = solve(prob, AutoTsit5(Rosenbrock23()))
end

main()
2

First of all, when it’s misbehaving use positive reinforcement instead of negative reinforcement. Give it a cookie for what it does well.

This is likely your issue right here because it’s a discontinuity. You might want to just make a continuous callback that hits that directly since otherwise the integrator can step over that assuming the simulation has continuous derivatives, which is violated at this point and makes its extrapolation inaccurate. That could lead to instability.