New to ModelingToolkit, problems trying to include derivatives on RHS of ODESystem

Specific Domains Modelling & Simulations
, ,
1

I’m trying to model a biological system of a population, N, and a trait, X.

The rate of change of the population, \frac{dN}{dt}, depends on some nonlinear function that is governed by some dynamic trait of the population, f(X).

The rate of change of X, \frac{dX}{dt}, depends on the derivative of \frac{dN}{dt} with respect to X: \frac{dX}{dt} = \frac{\partial}{\partial X}\frac{dN}{dt}.

I am trying to do this entirely in Julia, as my actual system is a bit more complex and I want to avoid introducing human error, meaning that I’m trying to take advantage of Symbolics differentiation. However, I run into problems when using the derivative on the right-hand side of the differential equations.

A demonstration:

 ## Setup demonstration
using ModelingToolkit, Symbolics, Plots

# I define a demo function here that is similar but simpler than what I actually
# need - a, b, and c are parameters, X a state variable
f(X, a, b, c) = a / (1 + (b / X) ^ c); 

# Usually I try it out with some test parameters to make sure I haven't done
# anything weird
p = (a = 1.0,
     b = 2.0,
     c = 0.5);

# And also with a test vector of X
Xvec = 0:0.01:2;

# Then I like to see how this function is behaving with X
plot(Xvec, f.(Xvec, p...), label = :none)

## Moving into ModelingToolkit stuff

# State variables and parameters
@parameters t a b c; # Including time, t
@variables N(t) X(t); # Two state variables to track, 
                      # the 2nd depends on the derivative w.r.t. X of the 1st

# I need a symbolic version of the earlier function for ModelingToolkit 
# to do its thing
F = f(X, a, b, c);

# Define derivatives
D = Differential(t); # time
Dx = Differential(X); # state variable X

# Rate of change of population
dNdt = N * F - 0.1;

# Defining the equations for the system
eqs = [D(N) ~ dNdt;
       D(X) ~ Dx(dNdt)]

# Putting it into ODE
@named sys = ODESystem(eqs); # this breaks
# returning: 
    # ArgumentError: Cannot find the parent of a / (1 + (b / X(t))^c).

Is there something specific that I’m doing wrong here or is this whole approach wrong? I am fairly new to this and haven’t found anything in the documentation that matches this exact problem.

Any help would be hugely appreciated !

2

Expand the derivatives so it does the symbolic differentiation:

@named sys = ODESystem(expand_derivatives.(eqs))

julia> equations(sys)
2-element Vector{Equation}:
 Differential(t)(N(t)) ~ (a*N(t)) / (1 + (b / X(t))^c)
 Differential(t)(X(t)) ~ (a*b*c*((b / X(t))^(c - 1))) / (((1 + (b / X(t))^c)^2)*(X(t)^2))
2 Likes
3

Amazing, thank you! I had tried expand_derivatives inside each line of eqs but not on the thing as a whole.

A follow-up question: is there a way to make this work for arrays of N and X? i.e., doing the same thing but now with the setting:

n = 2 # number of populations
@parameters t a[1:n] b[1:n] c[1:n]; # now parameters are vectors
@variables N(t)[1:n] X(t)[1:n]; # and so are the state variables
F = f.(X, a, b, c) # presumably need to broadcast function f to make it work
# and then what else?
4

It should just work? Do you have an MWE?

5

Sorry, was away for a few days.

The best attempt I have is this:

## NOW WITH ARRAYS OF STATE VARIABLES
using ModelingToolkit, Symbolics, DifferentialEquations, Plots;

# demo function again
f(X, a, b, c) = a / (1 + (b / X) ^ c); 

# State variables and parameters
@parameters t a[1:2] b[1:2] c[1:2]; # now as vectors
@variables N(t)[1:2] X(t)[1:2]; # now as vectors

# F function again
F = f.(X, a, b, c); # now needs broadcasting

# Define derivatives
D = Differential(t); # time
Dx = Differential(X); # state variable X

# Rate of change of population
dNdt = N .* F .- 0.1; # now needs broadcasting

# Defining the equations for the system
eqs = [D.(N) .~ dNdt; # now needs broadcasting
       D.(X) .~ Dx.(dNdt)];

# Define ODE system 
@named sys = ODESystem(expand_derivatives.(eqs)); 
 # ArgumentError: Cannot find the parent of (a[1]*(N(t))[1]) / (1 + (b[1] / (X(t))[1])^c[1]) - (X(t))[1].

I suppose I am doing some incorrect broadcasting?

6

The solution for this (thanks to @csimal) was to use Symbolics.scalarize.

# using scalarize inside the derivative of state variable X
Dx = Differential.(scalarize(X)) 
 
# creating dXdt by zipping together the "N"s and "Dx"s
dXdt = [d(x) for (d,x) in zip(Dx, Symbolics.scalarize(dNdt))]

# broadcasting the differentials on the result of scalarize
# everything gets turned into a vector
eqs = [D.(scalarize(N)) .~ scalarize(dNdt); # broadcasting "~"
       D.(scalarize(X)) .~ dXdt];

From here everything works smoothly!