How to create a birth/death system with ModelingToolkit.jl?

This gist of this question is to determine what the canonical workflow is for defining systems with changing numbers of components e.g. growth/death of biological organisms. This issue is discussed is the callback section of the DiffEq docs where the intended solution is to resize! the integrator according to your callback logic. This makes sense in the DE universe where you directly define julia functions which return the right hand side of your equation.

With MTK, it apparently makes its own rhs, so it wasn’t immediately clear to me how to define a model with a variable number of equations. More to the point, even if you can do such a thing, the resizing is happening at the level of the ODEProblem, so it’s not obvious what order your parameters are in anymore. One solution I cooked up is to decide a priori what the maximum number of things your system can have, and then multiply the rhs of all of your equations by an “activation” number which is in principle either 0 or 1 and is a parameter in the language of MTK. Then, when your callback logic triggers, you modify the appropriate activation variables to switch from 0 to 1 when something grows and vice versa when it dies. It still runs into the issue where you sort of need to know what order the parameters come in, but it seems like they’re in the order they’re defined? In general, maybe this isn’t a good idea since you’re pushing around possibly many extra equations which do nothing until they activate.

Here is a mwe for the activation method. I use Leaf as the thing that is growing/dying and a Tree as the glue between Leafs. To be clear, this doesn’t have anything to do with real trees and leaves. Of course, the ideal situation would be to be able to just add a new Leaf object directly in the callback logic and have everything Just Work™ but that might be tricky.

using ModelingToolkit, DifferentialEquations

struct LeafParameters{S<:Integer, T<:Real}

struct TreeParameters{T<:Real}

@variables t
ddt = Differential(t)

function Leaf(xy0::Vector{<:Real}, lp::LeafParameters; name::Symbol)
    ps = @parameters act = lp.activated xsp = lp.x_speed ysp = lp.y_speed
    @variables x(t) = xy0[1] y(t) = xy0[2]
    @variables F(t)[1:2]

    eqs = [
        ddt(x) ~ act*(xsp + F[1]), 
        ddt(y) ~ act*(ysp + F[2])

    return ODESystem(eqs, t, [x, y, F[1:2]...], ps; name)

function Tree(xy0::Matrix{<:Real}, lp::Vector{<:LeafParameters}, tp::TreeParameters; name::Symbol)
    # a tree is a composition of leaves
    N_leaves = size(xy0, 1)
    @named leaf 1:N_leaves i -> Leaf(xy0[i,:], lp[i])
    ps = @parameters tps = tp.stem_strength

    # in practice, these would be complicated interaction forces
    forces_x = [
        leaf[i].F[1] ~ tps for i = 1:N_leaves

    forces_y = [
        leaf[i].F[2] ~ tps for i = 1:N_leaves

    eqs = [forces_x ; forces_y]

    return compose(ODESystem(eqs, t, [], ps; name = name), leaf[1:N_leaves]...)   

n_leaves = 3

xy0 = [rand() for i = 1:n_leaves, j = 1:2]
lp = [
    LeafParameters(1.1, 2.1, 1),
    LeafParameters(1.2, 2.2, 0),
    LeafParameters(1.3, 2.3, 0)

tp = TreeParameters(3.0)

@named tree = Tree(xy0, lp, tp)
tree = structural_simplify(tree)

t_range = (0.0, 10.0)

prob = ODEProblem(

# prob.p = [3.0, 1.0, 1.1, 2.1, 0.0, 1.2, 2.2, 0.0, 1.3, 2.3]

function grow_leaf(t1::Real)
    condition(u, t, integrator) = t1 - t

    function affect!(integrator)
        u = integrator.u
        u[3] = 0.0 # x
        u[4] = 0.0 # y

        integrator.p[5] = 1.0 # 5 is the index of the second Leaf's activation

        return nothing

    return ContinuousCallback(condition, affect!)

@time sol = solve(prob, saveat = 1.0, callback = grow_leaf(3.0))

I don’t think resizing equations will work well with ModelingToolkit right now.

Should I read this as saying that the “activation parameter” solution is acceptable, given that resizing isn’t good in MTK? Or that I should avoid MTK altogether if I care about this specific kind of model?

Avoid MTK for these models for the time being.