I received this warning after doing a simulation:
> ┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
> └ @ OrdinaryDiffEq ~/.juliapro/JuliaPro_v1.3.1-3/packages/OrdinaryDiffEq/5igHh/src/solve.jl:455
> ┌ Warning: First function call produced NaNs. Exiting.
> └ @ OrdinaryDiffEq ~/.juliapro/JuliaPro_v1.3.1-3/packages/OrdinaryDiffEq/5igHh/src/initdt.jl:137
> ┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
> └ @ DiffEqBase ~/.juliapro/JuliaPro_v1.3.1-3/packages/DiffEqBase/XoVg5/src/integrator_interface.jl:323
Here is the model
# Model building
model = @model begin
@param begin
# Fixed effects
tIC50 ∈ RealDomain(lower = 0.0)
tIMAX ∈ RealDomain(lower = 0.0)
tMTT ∈ RealDomain(lower = 0.0)
tγ ∈ RealDomain(lower = 0.0)
tδ ∈ RealDomain(lower = 0.0)
tIT50 ∈ RealDomain(lower = 0.0)
tIMAT ∈ RealDomain(lower = 0.0)
tbase ∈ VectorDomain(2, lower = 0.0)
# occasions
# Microconstants
ICL ∈ ConstDomain(57.4)
IKA ∈ ConstDomain(0.763)
IBIO ∈ ConstDomain(0.317)
IQ3 ∈ ConstDomain(11.5)
IQ4 ∈ ConstDomain(39)
IV2 ∈ ConstDomain(36.2)
IV3 ∈ ConstDomain(561)
IV4 ∈ ConstDomain(54.4)
# Inter-individual variability
Ω ∈ PDiagDomain(9)
# Between-occasion variability.
Ω_occ ∈ RealDomain(lower = 0.0, init = 0.1)
Ω_occ_δ ∈ RealDomain(lower = 0.0, init = 0.1)
# Residuals models
σ_add ∈ RealDomain(lower = 0.0)
σ_prop ∈ RealDomain(lower = 0.0)
end
@random begin
η ~ MvNormal(Ω)
η_occ ~ MvNormal([Ω_occ, Ω_occ, Ω_occ])
η_occ_δ ~ MvNormal([Ω_occ_δ, Ω_occ_δ, Ω_occ_δ])
end
@covariates wt path cycl
@pre begin
IC50 = tIC50 * exp(η[1]+ η_occ[cycl]) #*occ1 + η_occ[2]*occ2 + η_occ[3]*occ3)
IMAX = tIMAX * exp(η[2])
MTT = tMTT * exp(η[3])
K = 4/MTT
γ = tγ * exp(η[4])
δ = tδ * exp(η[5]+ η_occ_δ[cycl]) #*occ1 + η_occ_δ[2]*occ2 + η_occ_δ[3]*occ3)
IT50 = tIT50 * exp(η[7])
IMAT = tIMAT * exp(η[8])
b0 = (tbase[1]* exp(η[6]) - (IMAT*(t/24))/(IT50+(t/24)))*(1 - path)
b1 = (tbase[2]* exp(η[9]))*path
base = b0 + b1
bioav = (prol = base,
tran1 = base,
tran2 = base,
tran3 = base,
circ = base
)
# Define microconstants
K67=IKA
K70=ICL/IV2
K78=IQ3/IV2
K87=IQ3/IV3
K79=IQ4/IV2
K97=IQ4/IV4
V7 =IV2/IBIO
#base = bas0 - (IMAT*(t/24))/(IT50+(t/24))
end
@vars begin
cp := Cent/V7
eff = (IMAX * cp)/(IC50+cp)
fbm = (base/circ)^γ
fbk = (base/circ)^δ
end
@dynamics begin
# PD
prol' = K * prol * fbm * (1-eff) - K*fbk*prol
tran1' = K * fbk* (prol - tran1)
tran2' = K * fbk* (tran1 - tran2)
tran3' = K * fbk* (tran2 - tran3)
circ' = K * fbk* tran3 - K*circ
# 3 compartments PK model
depot' = - K67 * depot
Cent' = K67 * depot + K87 * peri1 + K97*peri2 - Cent*(K78+K79+K70)
peri1' = K78 * Cent - K87 * peri1
peri2' = K79 * Cent - K97 * peri1
end
@derived begin
conc := @. Cent/V7
# combined error model
dv ~ @. Normal(conc, sqrt((conc*σ_prop)^2 + σ_add^2))
end
end
The command line after model building
# Simulations
# Parameter values are taken from the publication
param = (
# PD parameters
tIC50 = 0.0743,
tIMAX = 0.881,
tMTT = 95,
tγ = 0.498,
tδ = 0.177,
tbase = [195, 275],
tIT50 = 68.2,
tIMAT = 101,
# PK parameters
ICL =57.4,
IKA =0.763,
IBIO = 0.317,
IQ3 = 11.5,
IQ4 = 39,
IV2 = 36.2,
IV3 = 561,
IV4 = 54.4,
Ω = [0.663,0.0, 0.0871, 0.097, 0.228, 0.39,
0.0, 0.842, 0.0],
Ω_occ = 0.304,
Ω_occ_δ = 0.224,
σ_add = 0.135,
σ_prop = 0.19
)
# subjects
ev = DosageRegimen(60, time = 0, cmt =:depot)
using Random
Random.seed!(1234)
covar() = (wt = rand(55:110), cycl = rand([1, 2, 3]),
path = rand([0,1])) # 2 refers to 0
#population with covariates
step_size = 0.1
max_time = 48
pop = Population(map(i -> Subject(id=i,evs=ev,cvs=covar()), 1:125))
obs = simobs(model,pop,param,obstimes= 0:step_size:max_time)
Any thoughts?