ChrisRackauckas
commented
4 months ago Your packages are just old. On the latest version this seems fine:
using DifferentialEquations
using ModelingToolkit
schedule = 2
dosetimes_pk = [14, 28]
bolus_size = [10000.0,30000.0]
tspan = (0.0, dosetimes_pk[end]+10)
function pk_submodel_creator(;name, tinjection=[], bolus = 0.0)
@variables t
D = Differential(t)
ps = @parameters k_12=672 k_3=959 V_c=3410 V_p=3410
vars = @variables PK_c(t) PK_p(t)
bolus = length(bolus) < length(tinjection) ? repeat([bolus[1]], length(tinjection)) : bolus
injection = [[tinjection[j]] => [PK_c ~ PK_c + bolus[j]] for j in 1:length(tinjection)]
eqs = [
D(PK_c) ~ k_12 * PK_p/V_p - k_12 * PK_c/V_c - k_3 * PK_c/V_c,
D(PK_p) ~ k_12 * PK_c/V_c - k_12 * PK_p/V_p
]
ODESystem(eqs, t, vars, ps; name, discrete_events = injection)
end
function tcell_submodel_creator(;name)
@variables t
ps = @parameters σ, δ
vars = @variables E(t) ip(t)=0.0 [input = true] op(t)=0.0 [input = true]
D = Differential(t)
eqs = [
D(E) ~ σ - δ*E + ip - op,
]
ODESystem(eqs, t, vars, ps; name)
end
function tumourcell_submodel_creator(;name)
@variables t
ps = @parameters α
vars = @variables T(t) op(t)=0.0 [input = true] ip(t)=0.0 [input = true]
D = Differential(t)
eqs = [
D(T) ~ α * T - op + ip
]
ODESystem(eqs, t, vars, ps; name)
end
function extended_model(;name)
@named tcell = tcell_submodel_creator()
@named tumour = tumourcell_submodel_creator()
@named pk = pk_submodel_creator(; tinjection=dosetimes_pk, bolus=bolus_size)
@variables t
D = Differential(t)
ps = @parameters ρ, η, μ, c, ν, E_maxkill, EC_50_kill
mm_PK_kill = ((E_maxkill * pk.PK_c)/(EC_50_kill + pk.PK_c))
tumour_induced_tcell_proliferation = (ρ * tcell.E * (tumour.T)^c)/(η + (tumour.T)^c)
tumour_induced_tcell_death = μ * tcell.E * (tumour.T)^c
t_cell_induce_tumour_death = (ν + mm_PK_kill) * tcell.E * (tumour.T)^c
connections = [
tcell.ip ~ tumour_induced_tcell_proliferation,
tcell.op ~ tumour_induced_tcell_death,
tumour.op ~ t_cell_induce_tumour_death,
]
connected_blood_peripheral = compose(ODESystem(connections, t, [], ps; name), tcell, tumour, pk)
end
function create_ode_problem(model, tspan; kwargs...)
@variables t
u0 = [
model.tcell.E => 10e0,
model.tumour.T => 10e2,
model.pk.PK_c => 2500.0,
model.pk.PK_p => 0.0,
]
p = [
model.tumour.α => 0.12,
model.tcell.σ => 0.091,
model.tcell.δ => 0.00911,
model.ρ => 0.11,
model.η => 2.019e5,
model.μ => 3.422e-10,
model.c => 0.75,
model.ν => 1.101e-7,
model.E_maxkill => 4.801e-1,
model.EC_50_kill => 10000.0
]
return ODEProblem(model, u0, tspan, p; kwargs...)
end
@named model = extended_model();
model_simp = structural_simplify(model)
ode_prob = create_ode_problem(model_simp, tspan)
solve(ode_prob, Rosenbrock23())
solve(ode_prob)
solve(ode_prob; alg_hints=[:auto])
Describe the bug 🐞
The automated solver selection of
solveleads toDomainErrorwhile solving the system of ODEs and only explicitly setting a stiff solver such asRosenbrock23()leads to successfully solving the ODE problem.Expected behavior Speaking to @ChrisRackauckas on Slack, the solve with automated solver selection should not fail.
Minimal Reproducible Example 👇
Error & Stacktrace ⚠️ The last two
solvecalls lead to the following error:Environment:
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()