@wrholmes, welcome to the community! As @vijay said, combined pk/pd model computation is possible and (relatively) easy to accomplish. Consider the following 1CM for an IV bolus and a simple exponential response model (resp = E0 * (-kdrug * EXPOSURE)) where E0 is the baseline for some metric and kdrug is the slope of the drug effect. The EXPOSURE term can represent any PK metric computed by the model. I’ve provided 4 examples (resp1 to resp4) that I’ll discuss individually below. My focus will be on the content of the @derived block, which will change depending on which type of exposure you use to drive the PD response.
mdl = @model begin
@metadata begin
desc = "Combined PK/PD Model"
end
@param begin
#* PK
# Clearance, L/hr
tvcl ∈ RealDomain(init = 1.2)
# Central volume of distribution, L
tvvc ∈ RealDomain(init = 6.5)
# IIV
Ω ∈ PDiagDomain(init = [0.04,0.04])
# Proportional RUV
σpk ∈ RealDomain(init = 0.01)
#* PD
# baseline marker
tve0 ∈ RealDomain(init = 120)
# slope for drug effect
tvkdrug ∈ RealDomain(init = 0.2)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
cl = tvcl * exp(η[1])
vc = tvvc * exp(η[2])
e0 = tve0
kdrug = tvkdrug
end
@dynamics begin
central' = -cl/vc*central
# integrate of central compartment
int_central' = central/vc
end
@derived begin
#* PK
cp := @. central/vc
dv ~ @. Normal(cp, cp*σpk)
#* PD
# response driven by concentration
resp1 = @. e0 * exp(-kdrug * cp)
end
end
Response driven by concentration (resp1)
Take a look at the @dynamics, note that I’ve included a variable (int_central) which is used to integrate over central and provides a cumulative AUC. This variable will be used for all of the AUC related calculations that follow.
In the @derived block, I’ve included the standard cp (“IPRED”) and dv(“fitted observed values”), and I’ve used cp as my exposure metric to derive resp1. This application is fairly straightforward, and you can verify the result with a simple simobs call.
mysim = simobs(
mdl,
Subject(
events = DosageRegimen(50; addl = 2, ii = 24)
),
init_params(mdl),
zero_randeffs(mdl, init_params(mdl));
obstimes = 0:1:72
)
The init_params(mdl) function returns the parameter estimates provided in the @param block of the mdl definition. The zero_randeffs function sets the random effects (IIV) to 0. I’ve simulated hourly observations (0:1:72) for a dosage regimen of 50 “mass units” IV q24h x3.
From here, you’d just need to do a little post-processing to see the result.
Response driven by cumulative AUC (resp2)
If you’d like to use cumulative AUC, you’ll need to add it to the @derived block (the rest of the model is unchanged).
@derived begin
#* PK
cp := @. central/vc
dv ~ @. Normal(cp, cp*σpk)
# add cumulative auc
auc = @. int_central
#* PD
# response driven by cumulative AUC
resp2 = @. e0 * exp(-kdrug * auc)
end
The simobs call to test this is exactly the same as resp1.
Response driven by AUCtau (resp3)
(Note, for implementing resp3 and resp4 you need to call using ShiftedArrays: lead, lag in the REPL so you’ll have access to the lag function.)
I would guess that this response type is what you were referring to in your original post when you mentioned a “past-looking window.” The difference between this implementation and resp2 is the calculation of AUCtau using a combination of lag and coalesce (take a look at their individual docstrings for more detail). In short, we take the difference of int_central at any time t and the “lagged” vector we get by calling lag(int_central). The first lagged value is missing, thus I used coalesce to replace that missing value with 0 at time = 0. The broadcasting macro (@.) was omitted because you need to operate on the entire vector for lag to work.
@derived begin
#* PK
cp := @. central/vc
dv ~ @. Normal(cp, cp*σpk)
auctau = int_central .- coalesce.(lag(int_central), 0.0)
#* PD
# reponse driven by AUCτ (τ = q24h)
resp3 = @. e0 * exp(-kdrug * auctau)
end
The simobs call is also different for this implementation. Specifically, your obstimes are limited to 0 and the times corresponding to tau (i.e., multiples of your ii value). This is because the lag function would give nonsensical values for other times. For example, lag([0,1,...24])would give [missing, 0, 1, ...24] and AUCt1 - AUCt0 != AUCtau. The correct call would be:
mysim = simobs(
mdl,
Subject(
events = DosageRegimen(50; addl = 2, ii = 24)
),
init_params(mdl),
zero_randeffs(mdl, init_params(mdl));
obstimes = 0:24:72 #<- q24h
)
Response driven by average concentration (resp4)
The last implementation is similar to resp3 and has similar limitations. In this case, I’ve used “average concentration” by dividing the difference in int_central and lag(int_central) by an arbitrary value of 48 hours.
@derived begin
#* PK
cp := @. central/vc
dv ~ @. Normal(cp, cp*σpk)
cavg48 = (int_central .- coalesce.(lag(int_central), 0.0)) ./ 48
#* PD
# response driven by Cavg48
#resp4 = @. e0 * exp(-kdrug * cavg48)
end
For the simobs call, obstimes values need to be in 48-hour increments (or whatever period you choose). For this example, that would mean obstimes = [0,48], but if you were giving more doses, you could pass 0:48:X where X is some multiple of 48.
Hope that helps!
Best,
Haden