Hi all. Just trying to understand TMDD model. Here is the attached code that I am sharing. I am able to capture the graph of a typical TMDD profile, but I don’t replicate the plots in the book (Gabrielson PK27), I have tried my best to construct the model, I would highly appreciate someone to go through my model and identify the flaw…!
pk27 = @model begin
@param begin
tvcl ∈ RealDomain(lower=0)
tvkon ∈ RealDomain(lower=0)
tvkoff ∈ RealDomain(lower=0)
tvvc ∈ RealDomain(lower=0)
tvvp ∈ RealDomain(lower=0)
tvq ∈ RealDomain(lower=0)
tvkout ∈ RealDomain(lower=0)
tvkerl ∈ RealDomain(lower=0)
tvkin ∈ RealDomain(lower=0)
Ω ∈ PDiagDomain(9)
σ ∈ RealDomain(lower=0)
σ_prop ∈ RealDomain(lower=0)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
cl = tvcl * exp(η[1])
kon = tvkon * exp(η[2])
koff = tvkoff * exp(η[3])
vc = tvvc * exp(η[4])
vp = tvvp * exp(η[5])
q = tvq * exp(η[6])
kout = tvkout * exp(η[7])
kerl = tvkerl * exp(η[8])
kin = tvkin * exp(η[9])
ro = kin/kout
end
### I wanted to include ro = kin/kout in @init begin block, but when if do so,an error pops up saying ro not defined. So I had to define ro in the pre block! Why does that show up?
### Assuming my complex would be present in the central compartment, and it is the concentration of the complex, I use complex/vc in the equation !
@dynamics begin
Central' = - (cl/vc) * Central - (q/vc) * Central + (q/vp) * Peripheral - kon * (Central) * (ro) + koff * (Complex/vc)
Peripheral' = (q/vc) * Central - (q/vp) * Peripheral
Target' = kin - kout * (ro) - (kon) * Central * (ro) + koff * (Complex/vc)
Complex' = kon * (Central) * (ro) - (koff)* (Complex/vc) - kerl * (Complex/vc)
end
@derived begin
cp = @. Central/vc
dv ~ @. Normal(cp,sqrt(cp^2*σ))
target = @. Target
complex = @. Complex
end
end
Parameters
para = ( tvcl = 0.001,
tvkon = 0.096 ,
tvkoff = 0.001,
tvvc = 0.05,
tvvp = 0.100,
tvq = 0.003,
tvkout = 0.009,
tvro = 12,
tvkerl = 0.003,
tvkin = 0.11,
Ω = Diagonal([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]),
σ_prop = 0.00,
σ = 0.00)
evs1 = DosageRegimen(1.5,time=0,cmt=1)
S1 = Subject(id=1, evs=evs1)
evs2 = DosageRegimen(5,time=0,cmt=1)
S2 = Subject(id=2, evs=evs2)
evs3 = DosageRegimen(15,time=0,cmt=1)
S3 = Subject(id=3, evs=evs3)
evs4 = DosageRegimen(45,time=0,cmt=1)
S4 = Subject(id=4, evs=evs4)
population = Population([S1,S2,S3,S4])
obs1 = simobs(pk27,population,para,obstimes=0:0.1:500);
obs2 = simobs(pk27,population,para,obstimes=[0, 1, 10, 24, 72, 120, 168, 240, 360, 500])
obs3= simobs(pk27, population, para, obstimes=[0, 1, 5, 10, 10, 15, 20, 24, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 72, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 168, 170, 175, 180, 185, 190, 195, 200, 205, 210, 215, 220, 225, 230, 235, 240, 240, 245, 250, 255, 260, 265, 270, 275, 280, 285, 290, 295, 300, 305, 310, 315, 320, 325, 330, 335, 340, 345, 350, 355, 360, 360, 365, 370, 375, 380, 385, 390, 395, 400, 405, 410, 415, 420, 425, 430, 435, 440, 445, 450, 455, 460, 465, 470, 475, 480, 5, 490, 495, 500, 504])
df = DataFrame(obs1)
df1 = DataFrame(obs2)
df2 = DataFrame(obs3)
plot(obs1,yaxis=:log, yticks=[0.001,0.01,0.1,1,10,100,1000], ylims=(0.001,1000))
plot(obs2)
gr()
@df df plot(:time, :cp, yaxis=:log, yticks=[0.001,0.01,0.1,1,10,100,1000], ylims=(0.001,1000))
@df df plot!(:time, :target)
@df df plot!(:time, :complex)
@df df1 scatter!(:time, :target, color=[:blue])
@df df1 scatter!(:time, :complex, color=[:darkgreen])
@df df1 scatter!(:time, :cp, color=[:red])
