using Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
using CSV

Case Study III: Development of a population PKPD model
1 Introduction
This is the third tutorial in a series of case studies based on the tutorials found here . The third case study is about building a sequential PKPD model. It has IV infusion dosing, PK is governed by a simple one compartment model, and PD is an indirect response model (IDR) of histamine concentrations.
2 Data
The datasets are available from the PharmaDatasets
package. One for the PK model (CS3_IVINFEST
) and another for the PD model (CS3_IVINFPDEST
). First, we read the data for the PK model. We define an :evid
and a :cmt
column and set all event row values of :CONC
to missing
.
= CSV.read(dataset("pumas/event_data/CS3_IVINFEST", String), DataFrame; header = 4)
pkdata @rtransform!(pkdata, :evid = :AMT == 0 ? 0 : 1)
@rtransform!(pkdata, :cmt = :AMT == 0 ? missing : Symbol("Central"))
@rtransform!(pkdata, :CONC = :evid == 1 ? missing : :CONC)
Then, we map the DataFrame
to a population. We can omit specifying the cmt
and evid
keyword because we used the default value of lower case :cmt
and :evid
.
= read_pumas(
pk_population
pkdata;= :CID,
id = :TIME,
time = [:CONC],
observations = :AMT,
amt = :RATE,
rate )
Population
Subjects: 20
Observations: CONC
The data can be plotted using the observations_vs_time
plot.
observations_vs_time(pk_population; axis = (title = "PK data plot",))
To emphasize individual trajectories, the sim_plot
function also works on a population.
sim_plot(pk_population; axis = (title = "PK data plot",))
3 PK Model definition
The next step is to define the PK model. Since we have IV infusion we do not need a depot. The model does not contain any significant distribution phase. With just a Central
compartment, we can use the Central1
predefined model and its associated closed form solution. This is equivalent to ADVAN1
in NONMEM.
= @model begin
inf1cmt @param begin
∈ RealDomain(lower = 0.0)
θcl ∈ RealDomain(lower = 0.0)
θvc ∈ PDiagDomain(2)
Ω ∈ RealDomain(lower = 0.0)
σ_add ∈ RealDomain(lower = 0.0)
σ_prop end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= θcl * exp(η[1])
CL = θvc * exp(η[2])
Vc end
@dynamics Central1
@derived begin
:= @. Central / Vc
conc_model ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model * σ_prop)^2))
CONC end
end
PumasModel
Parameters: θcl, θvc, Ω, σ_add, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: CONC
Observed: CONC
To be able to fit the model we need to specify initial parameters.
= (
initial_est_inf1cmt = 1.0,
θcl = 1.0,
θvc = Diagonal([0.09, 0.09]),
Ω = sqrt(10.0),
σ_add = sqrt(0.01),
σ_prop )
(θcl = 1.0,
θvc = 1.0,
Ω = [0.09 0.0; 0.0 0.09],
σ_add = 3.1622776601683795,
σ_prop = 0.1,)
And then we can fit the model to data.
= fit(inf1cmt, pk_population, initial_est_inf1cmt, FOCE()) inf1cmt_results
[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 3.174194e+03 1.463122e+03 * time: 0.03469681739807129 1 1.700194e+03 3.899473e+02 * time: 1.7266778945922852 2 1.372117e+03 2.133637e+02 * time: 1.7342908382415771 3 1.149219e+03 9.540101e+01 * time: 1.742297887802124 4 1.063179e+03 5.172865e+01 * time: 1.7495667934417725 5 1.025366e+03 2.646221e+01 * time: 1.7566008567810059 6 1.011715e+03 1.666561e+01 * time: 1.7632169723510742 7 1.007591e+03 1.182557e+01 * time: 1.7694919109344482 8 1.006663e+03 9.681515e+00 * time: 1.7755868434906006 9 1.006394e+03 8.838980e+00 * time: 1.7815217971801758 10 1.006061e+03 8.916639e+00 * time: 1.78721284866333 11 1.005287e+03 8.898501e+00 * time: 1.7930629253387451 12 1.003673e+03 1.594601e+01 * time: 1.7990288734436035 13 1.000433e+03 3.264510e+01 * time: 1.8050928115844727 14 9.939620e+02 5.680776e+01 * time: 1.8112568855285645 15 9.920492e+02 8.189008e+01 * time: 1.8177649974822998 16 9.853470e+02 4.792287e+01 * time: 1.8244719505310059 17 9.821580e+02 4.099167e+01 * time: 1.8298218250274658 18 9.784445e+02 4.359738e+01 * time: 1.835981845855713 19 9.753844e+02 1.715016e+01 * time: 1.842167854309082 20 9.742523e+02 1.244736e+01 * time: 1.8471238613128662 21 9.734762e+02 1.086879e+01 * time: 1.8521909713745117 22 9.729031e+02 1.529972e+01 * time: 1.8571298122406006 23 9.714972e+02 2.051054e+01 * time: 1.8624320030212402 24 9.702763e+02 1.703802e+01 * time: 1.8676419258117676 25 9.693945e+02 1.242588e+01 * time: 1.8729209899902344 26 9.690627e+02 1.193813e+01 * time: 1.8779418468475342 27 9.687959e+02 1.049254e+01 * time: 1.8831830024719238 28 9.683343e+02 7.801000e+00 * time: 1.8882718086242676 29 9.667001e+02 1.096197e+01 * time: 1.8938570022583008 30 9.660129e+02 9.792112e+00 * time: 1.898535966873169 31 9.655534e+02 9.823330e+00 * time: 1.9033138751983643 32 9.651971e+02 1.059550e+01 * time: 1.9080917835235596 33 9.636185e+02 1.011640e+01 * time: 1.9135258197784424 34 9.621257e+02 7.727968e+00 * time: 1.9183719158172607 35 9.603739e+02 9.857407e+00 * time: 1.9234158992767334 36 9.603277e+02 1.508837e+01 * time: 1.928480863571167 37 9.598745e+02 3.570285e-01 * time: 1.933471918106079 38 9.598694e+02 2.341098e-01 * time: 1.9380009174346924 39 9.598690e+02 9.077958e-02 * time: 1.9424059391021729 40 9.598690e+02 3.615731e-02 * time: 1.9475657939910889 41 9.598690e+02 3.631911e-03 * time: 1.9527208805084229 42 9.598690e+02 3.934849e-04 * time: 1.9573948383331299
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 20
Observation records: Active Missing
CONC: 220 0
Total: 220 0
Number of parameters: Constant Optimized
0 6
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -959.86896
------------------
Estimate
------------------
θcl 0.024451
θvc 0.074289
Ω₁,₁ 0.072037
Ω₂,₂ 0.09384
σ_add 3.2726
σ_prop 0.10228
------------------
4 Model diagnostics
As usual, we use the inspect
function to calculate all diagnostics.
= inspect(inf1cmt_results) inf1cmt_insp
[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
FittedPumasModelInspection
Likelihood approximation used for weighted residuals: FOCE
These can be saved to a file by constructing a table representation of everything from predictions to weighted residuals and empirical bayes estimes and individual coefficients (POSTHOC in NONMEM).
= DataFrame(inf1cmt_insp)
df_inspect write("inspect_file.csv", df_inspect) CSV.
"inspect_file.csv"
Besides mean predictions, it is also simple to simulate from the estimated model using the empirical bayes estimates as the values for the random effects.
sim_plot(simobs(inf1cmt, pk_population, coef(inf1cmt_results)))
For the PD model we need individual CL
and Vc
values.
= unique(df_inspect[!, [:id, :time, :CL, :Vc]], :id)
icoef_dataframe rename!(icoef_dataframe, :CL => :CLi, :Vc => :Vci)
Row | id | time | CLi | Vci |
---|---|---|---|---|
String | Float64 | Float64? | Float64? | |
1 | 1 | 0.0 | 0.0154794 | 0.0879781 |
2 | 10 | 0.0 | 0.0308652 | 0.0566826 |
3 | 11 | 0.0 | 0.028236 | 0.0814447 |
4 | 12 | 0.0 | 0.0397456 | 0.0597957 |
5 | 13 | 0.0 | 0.021964 | 0.0764398 |
6 | 14 | 0.0 | 0.0233953 | 0.115324 |
7 | 15 | 0.0 | 0.0270355 | 0.0653006 |
8 | 16 | 0.0 | 0.0253244 | 0.0579478 |
9 | 17 | 0.0 | 0.0278282 | 0.144549 |
10 | 18 | 0.0 | 0.0400327 | 0.0813417 |
11 | 19 | 0.0 | 0.0183681 | 0.0723819 |
12 | 2 | 0.0 | 0.0191042 | 0.0920158 |
13 | 20 | 0.0 | 0.0276561 | 0.0533032 |
14 | 3 | 0.0 | 0.0156637 | 0.0753704 |
15 | 4 | 0.0 | 0.0260312 | 0.0691899 |
16 | 5 | 0.0 | 0.0284244 | 0.103891 |
17 | 6 | 0.0 | 0.0245165 | 0.094743 |
18 | 7 | 0.0 | 0.0153045 | 0.0360822 |
19 | 8 | 0.0 | 0.0261443 | 0.0750091 |
20 | 9 | 0.0 | 0.0248482 | 0.0555691 |
5 Getting the data
The data file consists of data obtained from 10 individuals who were treated with 500mg dose TID (three times a day, every eight hours) for five days. The dataset exists in the PharmaDatasets
package and we load it into memory as a DataFrame
. We specify that the first column should be a String
because we want to join the individual parameters from the PK step to the PD data, and the id
column of the DataFrame
from inspect
is a String
column. We use the first(input, number_of_elements)
function to show the first 10 rows of the DataFrame
.
= CSV.read(
pddata dataset("pumas/event_data/CS3_IVINFPDEST", String),
DataFrame;= 5,
header = Dict(1 => String),
types
)rename!(pddata, :CID => :id, :TIME => :time, :AMT => :amt, :CMT => :cmt)
@rtransform!(pddata, :evid = :amt == 0 ? 0 : 1)
@rtransform!(pddata, :HIST = :evid == 1 ? missing : :HIST)
first(pddata, 10)
Row | id | time | HIST | amt | cmt | RATE | MDV | CLI | VI | evid |
---|---|---|---|---|---|---|---|---|---|---|
String | Float64 | Float64? | Int64 | Int64 | Float64 | Int64 | Float64 | Float64 | Int64 | |
1 | 1 | 0.0 | missing | 100 | 1 | 16.7 | 1 | 15.585 | 90.812 | 1 |
2 | 1 | 0.0 | missing | 1 | 2 | 0.0 | 1 | 15.585 | 90.812 | 1 |
3 | 1 | 0.0 | 13.008 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
4 | 1 | 0.5 | 13.808 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
5 | 1 | 1.0 | 8.6859 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
6 | 1 | 3.0 | 6.2601 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
7 | 1 | 5.0 | 4.0602 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
8 | 1 | 6.0 | 4.4985 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
9 | 1 | 8.0 | 3.2736 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
10 | 1 | 12.0 | 2.6026 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 |
The file contains 11 columns:
= outerjoin(pddata, icoef_dataframe; on = [:id, :time])
pd_dataframe sort!(pd_dataframe, [:id, :time])
Row | id | time | HIST | amt | cmt | RATE | MDV | CLI | VI | evid | CLi | Vci |
---|---|---|---|---|---|---|---|---|---|---|---|---|
String | Float64 | Float64? | Int64? | Int64? | Float64? | Int64? | Float64? | Float64? | Int64? | Float64? | Float64? | |
1 | 1 | 0.0 | missing | 100 | 1 | 16.7 | 1 | 15.585 | 90.812 | 1 | 0.0154794 | 0.0879781 |
2 | 1 | 0.0 | missing | 1 | 2 | 0.0 | 1 | 15.585 | 90.812 | 1 | 0.0154794 | 0.0879781 |
3 | 1 | 0.0 | 13.008 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | 0.0154794 | 0.0879781 |
4 | 1 | 0.5 | 13.808 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
5 | 1 | 1.0 | 8.6859 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
6 | 1 | 3.0 | 6.2601 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
7 | 1 | 5.0 | 4.0602 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
8 | 1 | 6.0 | 4.4985 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
9 | 1 | 8.0 | 3.2736 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
10 | 1 | 12.0 | 2.6026 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
11 | 1 | 15.0 | 2.3422 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
12 | 1 | 18.0 | 3.8008 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
13 | 1 | 24.0 | 7.1397 | 0 | 2 | 0.0 | 0 | 15.585 | 90.812 | 0 | missing | missing |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
269 | 9 | 0.0 | 34.608 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | 0.0248482 | 0.0555691 |
270 | 9 | 0.5 | 29.57 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
271 | 9 | 1.0 | 27.68 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
272 | 9 | 3.0 | 18.397 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
273 | 9 | 5.0 | 14.908 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
274 | 9 | 6.0 | 11.669 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
275 | 9 | 8.0 | 8.5475 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
276 | 9 | 12.0 | 16.895 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
277 | 9 | 15.0 | 21.326 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
278 | 9 | 18.0 | 26.131 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
279 | 9 | 24.0 | 26.278 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
280 | 9 | 30.0 | 39.643 | 0 | 2 | 0.0 | 0 | 24.809 | 56.453 | 0 | missing | missing |
6 Converting the DataFrame to a collection of Subjects
= read_pumas(pd_dataframe; observations = [:HIST], covariates = [:CLi, :Vci]) population
Population
Subjects: 20
Covariates: CLi, Vci
Observations: HIST
7 IDR model
= @model begin
irm1 @metadata begin
= "POPULATION PK-PD MODELING"
desc = u"hr" # hour
timeu end
@param begin
∈ RealDomain(lower = 0)
tvkin ∈ RealDomain(lower = 0)
tvkout ∈ RealDomain(lower = 0)
tvic50 ∈ RealDomain(lower = 0)
tvimax ∈ PDiagDomain(3)
Ω ∈ RealDomain(lower = 0)
σ_add_pd ∈ RealDomain(lower = 0)
σ_prop_pd end
@random begin
~ MvNormal(Ω)
η end
@covariates CLi Vci
@pre begin
= tvkin * exp(η[1])
kin = tvkout * exp(η[2])
kout = kin / kout
bsl = tvic50 * exp(η[3])
ic50 = tvimax
imax = CLi
CL = Vci
Vc end
@init begin
= bsl
Response end
@dynamics begin
' = -CL / Vc * Central
Central' =
Response* (1 - imax * (Central / Vc) / (ic50 + Central / Vc)) - kout * Response
kin end
@derived begin
~ @. Normal(Response, sqrt(σ_add_pd^2 + (Response * σ_prop_pd)^2))
HIST end
end
PumasModel
Parameters: tvkin, tvkout, tvic50, tvimax, Ω, σ_add_pd, σ_prop_pd
Random effects: η
Covariates: CLi, Vci
Dynamical system variables: Central, Response
Dynamical system type: Nonlinear ODE
Derived: HIST
Observed: HIST
= (
init_θ = 5.4,
tvkin = 0.3,
tvkout = 3.9,
tvic50 = 1.0,
tvimax = Diagonal([0.2, 0.2, 0.2]),
Ω = 0.05,
σ_add_pd = 0.05,
σ_prop_pd )
(tvkin = 5.4,
tvkout = 0.3,
tvic50 = 3.9,
tvimax = 1.0,
Ω = [0.2 0.0 0.0; 0.0 0.2 0.0; 0.0 0.0 0.2],
σ_add_pd = 0.05,
σ_prop_pd = 0.05,)
= fit(
irm1_results
irm1,
population,
init_θ,FOCE();
Pumas.= (tvimax = 1.0,),
constantcoef = (show_every = 10,),
optim_options )
[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 1.625453e+03 2.073489e+03 * time: 1.811981201171875e-5 10 5.899604e+02 6.134653e+00 * time: 1.6497101783752441 20 5.824375e+02 2.641109e+01 * time: 2.58089017868042 30 5.801204e+02 1.223989e+00 * time: 3.4438860416412354 40 5.798576e+02 5.986551e-01 * time: 4.279842138290405 50 5.668948e+02 6.290186e+00 * time: 5.257589101791382 60 5.591291e+02 5.129868e+00 * time: 6.195143222808838 70 5.559638e+02 3.219738e-04 * time: 7.0542311668396
FittedPumasModel
Dynamical system type: Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)
Number of subjects: 20
Observation records: Active Missing
HIST: 240 0
Total: 240 0
Number of parameters: Constant Optimized
1 8
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -555.96385
------------------------
Estimate
------------------------
tvkin 5.5533
tvkout 0.27864
tvic50 30.427
† tvimax 1.0
Ω₁,₁ 0.18348
Ω₂,₂ 0.067395
Ω₃,₃ 0.20683
σ_add_pd 1.1094
σ_prop_pd 0.1038
------------------------
† indicates constant parameters
= inspect(irm1_results)
irm1_insp goodness_of_fit(irm1_insp)
[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
8 Conclusion
In this tutorial, we saw how to build on topics we learnt in the previous two case studies to build a sequential PKPD model. We built two different models and saw how to forward the results from the PK model to the PD model. This concludes the third of the three case studies.