Case Study III: Development of a population PKPD model

Author

Patrick Kofod Mogensen 

using Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
using CSV

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 (inf_sd_4_pk) and another for the PD model (inf_sd_4_idr). 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.

pkdata = dataset("inf_sd_4_pk")
@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.

pk_population = read_pumas(
    pkdata;
    id = :CID,
    time = :TIME,
    observations = [:CONC],
    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.

inf1cmt = @model begin
    @param begin
        θcl ∈ RealDomain(lower = 0.0)
        θvc ∈ RealDomain(lower = 0.0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(lower = 0.0)
        σ_prop ∈ RealDomain(lower = 0.0)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @pre begin
        CL = θcl * exp(η[1])
        Vc = θvc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
        conc_model := @. Central / Vc
        CONC ~ @. CombinedNormal(conc_model, σ_add, σ_prop)
    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 = (
    θcl = 1.0,
    θvc = 1.0,
    Ω = Diagonal([0.09, 0.09]),
    σ_add = sqrt(10.0),
    σ_prop = sqrt(0.01),
)
(θ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.

inf1cmt_results = fit(inf1cmt, pk_population, initial_est_inf1cmt, FOCE())
[ 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.033746957778930664
     1     1.700194e+03     3.899473e+02
 * time: 2.600399971008301
     2     1.372117e+03     2.133637e+02
 * time: 2.6078529357910156
     3     1.149219e+03     9.540101e+01
 * time: 2.6158039569854736
     4     1.063179e+03     5.172865e+01
 * time: 2.6232779026031494
     5     1.025366e+03     2.646221e+01
 * time: 2.6303908824920654
     6     1.011715e+03     1.666561e+01
 * time: 2.637495994567871
     7     1.007591e+03     1.182557e+01
 * time: 2.644467830657959
     8     1.006663e+03     9.681515e+00
 * time: 2.6512458324432373
     9     1.006394e+03     8.838980e+00
 * time: 2.6577727794647217
    10     1.006061e+03     8.916639e+00
 * time: 2.6645147800445557
    11     1.005287e+03     8.898501e+00
 * time: 2.670984983444214
    12     1.003673e+03     1.594601e+01
 * time: 2.677963972091675
    13     1.000433e+03     3.264510e+01
 * time: 2.6851229667663574
    14     9.939620e+02     5.680776e+01
 * time: 2.6923367977142334
    15     9.920492e+02     8.189008e+01
 * time: 2.699016809463501
    16     9.853470e+02     4.792287e+01
 * time: 2.705772876739502
    17     9.821580e+02     4.099167e+01
 * time: 2.7120108604431152
    18     9.784445e+02     4.359738e+01
 * time: 2.718787908554077
    19     9.753844e+02     1.715016e+01
 * time: 2.7253668308258057
    20     9.742523e+02     1.244736e+01
 * time: 2.731024980545044
    21     9.734762e+02     1.086879e+01
 * time: 2.7368969917297363
    22     9.729031e+02     1.529972e+01
 * time: 2.7420568466186523
    23     9.714972e+02     2.051054e+01
 * time: 2.747314929962158
    24     9.702763e+02     1.703802e+01
 * time: 2.7537078857421875
    25     9.693945e+02     1.242588e+01
 * time: 2.76005482673645
    26     9.690627e+02     1.193813e+01
 * time: 2.7663419246673584
    27     9.687959e+02     1.049254e+01
 * time: 2.7726149559020996
    28     9.683343e+02     7.801000e+00
 * time: 2.778930902481079
    29     9.667001e+02     1.096197e+01
 * time: 2.7853968143463135
    30     9.660129e+02     9.792112e+00
 * time: 2.883004903793335
    31     9.655534e+02     9.823330e+00
 * time: 2.888199806213379
    32     9.651971e+02     1.059550e+01
 * time: 2.893094778060913
    33     9.636185e+02     1.011640e+01
 * time: 2.898237943649292
    34     9.621257e+02     7.727968e+00
 * time: 2.903275966644287
    35     9.603739e+02     9.857406e+00
 * time: 2.908641815185547
    36     9.603277e+02     1.508837e+01
 * time: 2.914078950881958
    37     9.598745e+02     3.570285e-01
 * time: 2.919501781463623
    38     9.598694e+02     2.341098e-01
 * time: 2.9245519638061523
    39     9.598690e+02     9.077958e-02
 * time: 2.9293129444122314
    40     9.598690e+02     3.615731e-02
 * time: 2.9339327812194824
    41     9.598690e+02     3.631910e-03
 * time: 2.938560962677002
    42     9.598690e+02     3.934849e-04
 * time: 2.942922830581665

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.

inf1cmt_insp = inspect(inf1cmt_results)
[ 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).

df_inspect = DataFrame(inf1cmt_insp)
CSV.write("inspect_file.csv", df_inspect)
"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.

icoef_dataframe = unique(df_inspect[!, [:id, :time, :CL, :Vc]], :id)
rename!(icoef_dataframe, :CL => :CLi, :Vc => :Vci)
20×4 DataFrame
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 convert the subject identification column to 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.

pddata = dataset("inf_sd_4_idr")
pddata.CID = string.(pddata.CID)
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)
10×10 DataFrame
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:

pd_dataframe = outerjoin(pddata, icoef_dataframe; on = [:id, :time])
sort!(pd_dataframe, [:id, :time])
280×12 DataFrame
255 rows omitted
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

population = read_pumas(pd_dataframe; observations = [:HIST], covariates = [:CLi, :Vci])
Population
  Subjects: 20
  Covariates: CLi, Vci
  Observations: HIST

7 IDR model

irm1 = @model begin
    @metadata begin
        desc = "POPULATION PK-PD MODELING"
        timeu = u"hr" # hour
    end
    @param begin
        tvkin ∈ RealDomain(lower = 0)
        tvkout ∈ RealDomain(lower = 0)
        tvic50 ∈ RealDomain(lower = 0)
        tvimax ∈ RealDomain(lower = 0)
        Ω ∈ PDiagDomain(3)
        σ_add_pd ∈ RealDomain(lower = 0)
        σ_prop_pd ∈ RealDomain(lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates CLi Vci

    @pre begin
        kin = tvkin * exp(η[1])
        kout = tvkout * exp(η[2])
        bsl = kin / kout
        ic50 = tvic50 * exp(η[3])
        imax = tvimax
        CL = CLi
        Vc = Vci
    end

    @init begin
        Response = bsl
    end

    @dynamics begin
        Central' = -CL / Vc * Central
        Response' =
            kin * (1 - imax * (Central / Vc) / (ic50 + Central / Vc)) - kout * Response
    end

    @derived begin
        HIST ~ @. CombinedNormal(Response, σ_add_pd, σ_prop_pd)
    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_θ = (
    tvkin = 5.4,
    tvkout = 0.3,
    tvic50 = 3.9,
    tvimax = 1.0,
    Ω = Diagonal([0.2, 0.2, 0.2]),
    σ_add_pd = 0.05,
    σ_prop_pd = 0.05,
)
(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)
irm1_results = fit(
    irm1,
    population,
    init_θ,
    Pumas.FOCE();
    constantcoef = (:tvimax,),
    optim_options = (show_every = 10,),
)
[ 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.9788742065429688e-5
    10     5.899604e+02     6.134653e+00
 * time: 1.8663599491119385
    20     5.824375e+02     2.641119e+01
 * time: 2.4166557788848877
    30     5.801204e+02     1.223993e+00
 * time: 3.001513957977295
    40     5.798576e+02     5.986549e-01
 * time: 3.559739828109741
    50     5.668946e+02     6.290059e+00
 * time: 4.209833860397339
    60     5.591269e+02     5.067801e+00
 * time: 4.831350803375244
    70     5.559638e+02     2.302270e-03
 * time: 5.461353778839111
    80     5.559638e+02     1.418130e-03
 * time: 6.2952799797058105

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:              NoObjectiveChange
Log-likelihood value:                   -555.96384

------------------------
              Estimate
------------------------
  tvkin        5.5533
  tvkout       0.27864
  tvic50      30.427
† tvimax       1.0
  Ω₁,₁         0.18348
  Ω₂,₂         0.067394
  Ω₃,₃         0.20683
  σ_add_pd     1.1094
  σ_prop_pd    0.1038
------------------------
† indicates constant parameters
irm1_insp = inspect(irm1_results)
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.