Sequential and Simultaneous PKPD Models

Author

Vijay Ivaturi

In this tutorial, the approach to both sequential and simultaneous PK/PD modeling in Pumas using the Warfarin example will be demonstrated. The tutorial will:

Define and fit the previously established Warfarin PK model.
Demonstrate how to set up and fit a PD (PCA endpoint) model sequentially using the empirical Bayes estimates (EBEs) from the PK model.
Demonstrate how to set up and fit a simultaneous PK/PD model.
Discuss advantages/disadvantages of each approach.

1 Introduction

When working on pharmacokinetic/pharmacodynamic (PK/PD) problems, there are two broad choices:

Sequential modeling: Fit the PK model first and generate individual PK parameters. Then pass these individual parameters (as if they were data) into the PD model. This is sometimes called the “two-stage” or “sequential” approach.
Simultaneous (or joint) modeling: Fit the PK and PD data within one model, estimating PK and PD parameters at the same time.

1.1 Advantages and Disadvantages

Sequential approach

Advantages:
- Simpler, conceptually.
- Faster,
  - because the PK fit is separated from the PD fit.
  - if PD is non-linear and PK is linear at which point one can use closed form or matrix exponential and investigate the PK form much faster than if everything was fitted together.
- Easier to troubleshoot if the PD portion misbehaves (the PK portion is already finalized).
Disadvantages:
- The uncertainty from the PK step is not propagated into the PD parameter estimates.
- The PK model parameters are “frozen” once the PD step starts, which might overlook any interplay or correlation between PK and PD parameters.

Simultaneous approach

Advantages:
- A single model handles PK and PD data together, ensuring that the correlation and covariance between PK parameters and PD parameters are accounted for.
- Proper propagation of parameter uncertainty from PK into PD because the entire system is fit at once.
Disadvantages:
- Model can become more complex and harder to fit.
- Potentially heavier computational burden.
- More complicated model diagnostics.

Note

Readers familiar with the PK model building from the previous tutorial can choose to skip the PK model fitting section and proceed to the PKPD modeling section (Section 3).

2 The Warfarin PK Model

Let us reuse the previously established Warfarin PK model.

warfarin_pk_model = @model begin
    @metadata begin
        desc = "Warfarin 1-compartment PK model (PD removed)"
        timeu = u"hr"
    end
    @param begin
        # PK parameters
        """
        Clearance (L/hr)
        """
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central volume (L)
        """
        pop_Vc ∈ RealDomain(lower = 0.0, init = 8.11)
        """
        Absorption lag time (hr)
        """
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
        """
        Lag time (hr)
        """
        pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)
        # Inter-individual variability
        """
          - ΩCL: Clearance
          - ΩVc: Central volume
          - Ωtabs: Absorption lag time
        """
        pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])
        # Residual variability
        """
        σ_prop: Proportional error
        """
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        """
        σ_add: Additive error
        """
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
    end
    @random begin
        pk_η ~ MvNormal(pk_Ω)    # mean = 0, covariance = pk_Ω
    end
    @covariates begin
        """
        FSZCL: Clearance scaling factor
        """
        FSZCL
        """
        FSZV: Volume scaling factor
        """
        FSZV
    end
    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
    end
    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end
    @vars begin
        cp := Central / Vc
    end
    @dynamics Depots1Central1

    @derived begin
        """
        Concentration (ng/mL)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
    end
end

PumasModel
  Parameters: pop_CL, pop_Vc, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add
  Random effects: pk_η
  Covariates: FSZCL, FSZV
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: conc
  Observed: conc

# Source the example dataset
warfarin_data = dataset("paganz2024/warfarin_long")
# Some subjects have duplicate time points for DVID = 1
# For this dataset, the triple (ID, TIME, DVID) should define
# a row uniquely, but
nrow(warfarin_data)
nrow(unique(warfarin_data, ["ID", "TIME", "DVID"]))
# We can identify the problematic rows by grouping on the index variables
@chain warfarin_data begin
    @groupby :ID :TIME :DVID
    @transform :tmp = length(:ID)
    @rsubset :tmp > 1
end
# It is important to understand the reason for the duplicate values.
# Sometimes the duplication is caused by recording errors, sometimes
# it is a data processing error, e.g. when joining tables, or it can
# be genuine records, e.g. when samples have been analyzed in multiple
# labs. The next step depends on which of the causes are behind the
# duplications.
#
# In this case, we will assume that both values are informative and
# we will therefore just adjust the time stamp a bit for the second
# observation.
warfarin_data_processed = @chain warfarin_data begin
    @groupby :ID :TIME :DVID
    @transform :tmp = 1:length(:ID)
    @rtransform :TIME = :tmp == 1 ? :TIME : :TIME + 1e-6
    @select Not(:tmp)
end
# Transform the data in a single chain of operations
warfarin_data_wide = @chain warfarin_data_processed begin
    @rsubset !contains(:ID, "#")
    @rtransform begin
        # Scaling factors
        :FSZV = :WEIGHT / 70            # volume scaling
        :FSZCL = (:WEIGHT / 70)^0.75     # clearance scaling (allometric)
        # Column name for the DV
        :DVNAME = "DV$(:DVID)"
        # Dosing indicator columns
        :CMT = ismissing(:AMOUNT) ? missing : 1
        :EVID = ismissing(:AMOUNT) ? 0 : 1
    end
    unstack(Not([:DVID, :DVNAME, :DV]), :DVNAME, :DV)
    rename!(:DV1 => :conc, :DV2 => :pca)
end

pk_pop = read_pumas(
    warfarin_data_wide;
    id = :ID,
    time = :TIME,
    amt = :AMOUNT,
    cmt = :CMT,
    evid = :EVID,
    covariates = [:SEX, :WEIGHT, :FSZV, :FSZCL],
    observations = [:conc],
)

pkpd_pop = read_pumas(
    warfarin_data_wide;
    id = :ID,
    time = :TIME,
    amt = :AMOUNT,
    cmt = :CMT,
    evid = :EVID,
    covariates = [:SEX, :WEIGHT, :FSZV, :FSZCL],
    observations = [:conc, :pca],
)

Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL
  Observations: conc, pca

# A named tuple of parameter values
param_vals = (
    pop_CL = 0.134,
    pop_Vc = 8.11,
    pop_tabs = 0.523,
    pop_lag = 0.1,
    pk_Ω = Diagonal([0.01, 0.01, 0.01]),
    σ_prop = 0.00752,
    σ_add = 0.0661,
)


pk_fit = fit(warfarin_pk_model, pk_pop, param_vals, FOCE(););

FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:                 Closed form

Log-likelihood value:                   -319.93444
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0              9
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

-----------------------
             Estimate
-----------------------
pop_CL        0.13065
pop_Vc        8.0713
pop_tabs      0.37143
pop_lag       0.92611
pk_Ω₁,₁       0.050903
pk_Ω₂,₂       0.019133
pk_Ω₃,₃       1.1401
σ_prop        0.090622
σ_add         0.33356
-----------------------

3 Sequential PK/PD Modeling

In a sequential approach, the PD endpoint (e.g., PCA or clotting time, or INR) is modeled as a function of the Warfarin concentration using the indirect response model.

\frac{dTurnover}{dt} = R_{in} \cdot (1 + E_{max} \cdot \frac{C_{p}}{EC_{50} + C_{p}}) - k_{out} \cdot Turnover

The standard workflow is:

Fit the PK model (as above).
Extract the individual empirical Bayes estimates (EBEs) for each subject.
Create a new dataset for the PD model that includes these individual coefficients as inputs.
Read the new dataset into a Pumas Population object.
Fit the PD model to the PD population.

3.1 Extracting Individual Empirical Bayes Estimates

There are three ways to extract the individual EBEs:

Extract the individual coefficients (e.g., CL, Vc, tabs, etc.) from the fit and merge them with the original data using data wrangling tools.
Use the Subject constructor that augments the original PKPD data with the individual coefficients.

These methods require the full PKPD model to be defined and use the individual coefficients from the PK fit to predict the PK concentrations that get directly passed to the PD model. Below, the PKPD model is defined as warfarin_pkpd_model.

warfarin_pkpd_model = @model begin

    @metadata begin
        desc = "Warfarin PK/PD model"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Clearance (L/h/70kg)
        """
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central Volume L/70kg
        """
        pop_V ∈ RealDomain(lower = 0.0, init = 8.11)
        """
        Absorption time (h)
        """
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
        """
        Lag time (h)
        """
        pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)

        # PD parameters
        """
        Baseline
        """
        pop_e0 ∈ RealDomain(lower = 0.0, init = 100.0)
        """
        Emax
        """
        pop_emax ∈ RealDomain(init = -1.0)
        """
        EC50
        """
        pop_c50 ∈ RealDomain(lower = 0.0, init = 1.0)
        """
        Turnover
        """
        pop_tover ∈ RealDomain(lower = 0.0, init = 14.0)
        # Inter-individual variability
        """
          - ΩCL
          - ΩVc
          - ΩTabs
        """
        pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])

        """
          - Ωe0
          - Ωemax
          - Ωec50
          - Ωturn
        """
        pd_Ω ∈ PDiagDomain([0.01, 0.01, 0.01, 0.01])
        # Residual variability
        """
        Proportional residual error for drug concentration
        """
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        """
        Additive residual error for drug concentration (mg/L)
        """
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
        """
        Additive error for PCA
        """
        σ_fx ∈ RealDomain(lower = 0.0, init = 0.01)
    end

    @random begin
        # mean = 0, covariance = pk_Ω
        pk_η ~ MvNormal(pk_Ω)
        # mean = 0, covariance = pd_Ω
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates begin
        """
        Scaling factor for Volume
        """
        FSZV
        """
        Scaling factor for Clearance
        """
        FSZCL
    end

    @pre begin
        # PK
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_V * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs

        # PD
        e0 = pop_e0 * exp(pd_η[1])
        emax = pop_emax * exp(pd_η[2])
        c50 = pop_c50 * exp(pd_η[3])
        tover = pop_tover * exp(pd_η[4])
        kout = log(2) / tover
        rin = e0 * kout
    end

    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end

    @init begin
        Turnover = e0
    end

    @vars begin
        cp := Central / Vc
        ratein := Ka * Depot
        pd := 1 + emax * cp / (c50 + cp)
    end

    @dynamics begin
        Depot' = -ratein
        Central' = ratein - CL * cp
        Turnover' = rin * pd - kout * Turnover
    end

    @derived begin
        """
        Warfarin Concentration (mg/L)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
        """
        PCA
        """
        pca ~ @. Normal(Turnover, σ_fx)
    end
end

PumasModel
  Parameters: pop_CL, pop_V, pop_tabs, pop_lag, pop_e0, pop_emax, pop_c50, pop_tover, pk_Ω, pd_Ω, σ_prop, σ_add, σ_fx
  Random effects: pk_η, pd_η
  Covariates: FSZV, FSZCL
  Dynamical system variables: Depot, Central, Turnover
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca
  Observed: conc, pca

In this model above, the @pre block contains the expressions for the PK and PD parameters. The intention in a sequential approach is to let the individual coefficients of the PK model be fixed and known, and only estimate the PD parameters.

To facilitate this, the individual coefficients from the PK fit are introduced as @covariates in the @pre block. Hence, blocks representing the PK parameters such as @param, @random, @covariates, and the @pre blocks will be modified from the above specification to include the individual coefficients as seen below.

@param begin
    # Notice how all PK related parameters are now removed from the @param block

    # PD parameters
    """
    Baseline
    """
    pop_e0 ∈ RealDomain(lower = 0.0, init = 100.0)
    """
    Emax
    """
    pop_emax ∈ RealDomain(init = -1.0)
    """
    EC50
    """
    pop_c50 ∈ RealDomain(lower = 0.0, init = 1.0)
    """
    Turnover
    """
    pop_tover ∈ RealDomain(lower = 0.0, init = 14.0)
    # PD Inter-individual variability
    """
      - Ωe0
      - Ωemax
      - Ωec50
      - Ωturn
    """
    pd_Ω ∈ PDiagDomain([0.01, 0.01, 0.01, 0.01])
    # Residual variability
    """
    Proportional residual error for drug concentration
    """
    σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
    """
    Additive residual error for drug concentration (mg/L)
    """
    σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
    """
    Additive error for PCA
    """
    σ_fx ∈ RealDomain(lower = 0.0, init = 0.01)
end

@random begin
    # Notice how the random effects are now defined in terms of the PD parameters only
    # mean = 0, covariance = pd_Ω
    pd_η ~ MvNormal(pd_Ω)
end

@covariates begin
    # Notice how the individual coefficients are now defined in terms of the covariates
    """
    Individual Clearance
    """
    CL
    """
    Individual Volume
    """
    Vc
    """
    Individual Absorption time
    """
    tabs
    """
    Individual Lag time
    """
    lags_Depot
    """
    Scaling factor for Volume
    """
    FSZV
    """
    Scaling factor for Clearance
    """
    FSZCL
end

@pre begin
    # PK
    Ka = log(2) / tabs

    # PD
    e0 = pop_e0 * exp(pd_η[1])
    emax = pop_emax * exp(pd_η[2])
    c50 = pop_c50 * exp(pd_η[3])
    tover = pop_tover * exp(pd_η[4])
    kout = log(2) / tover
    rin = e0 * kout
end

@dosecontrol begin
    lags = (Depot = lags_Depot,)
end

To bring the individual coefficients into the model, we need to extract the individual coefficients from the PK fit and then introduce them into the original PKPD dataset or in the PKPD population. Below is an example of how to do this using both methods.

Step 1: Extract the individual coefficients from the PK fit.

# Perform inspection on FittedPumasModel
warf_pk_insp = inspect(pk_fit);
# Convert to a DataFrame
df_inspect = DataFrame(warf_pk_insp);
# Obtain the individual PK parameters
icoef_dataframe = unique(df_inspect[!, [:id, :time, :CL, :Vc, :tabs, :lags_Depot]], :id)
rename!(icoef_dataframe, :id => :ID, :time => :TIME);
first(icoef_dataframe, 5)

[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.

5×6 DataFrame

Row	ID	TIME	CL	Vc	tabs	lags_Depot
	String	Float64	Float64?	Float64?	Float64?	Float64?
1	1	0.0	0.112039	7.67823	0.379683	0.926108
2	11	0.0	0.131106	7.7971	0.313889	0.926108
3	12	0.0	0.134624	6.83003	2.6345	0.926108
4	13	0.0	0.174219	8.78166	0.829443	0.926108
5	14	0.0	0.152349	6.15482	0.209459	0.926108

Step 2: Merge the individual coefficients with the original PKPD dataset.

# Merge the individual PK parameters with the original dataset by ID and TIME
pd_dataframe = outerjoin(warfarin_data_wide, icoef_dataframe; on = [:ID, :TIME]);
# Arrange the DataFrame by ID and TIME
sort!(pd_dataframe, [:ID, :TIME]);
# Return the first 5 rows
first(pd_dataframe, 5)

5×17 DataFrame

Row	ID	TIME	WEIGHT	AGE	SEX	AMOUNT	FSZV	FSZCL	CMT	EVID	DV0	pca	conc	CL	Vc	tabs	lags_Depot
	String	Float64	Float64?	Int64?	Int64?	Float64?	Float64?	Float64?	Int64?	Int64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?
1	1	0.0	66.7	50	1	100.0	0.952857	0.96443	1	1	missing	missing	missing	0.112039	7.67823	0.379683	0.926108
2	1	0.0	66.7	50	1	missing	0.952857	0.96443	missing	0	missing	100.0	missing	0.112039	7.67823	0.379683	0.926108
3	1	24.0	66.7	50	1	missing	0.952857	0.96443	missing	0	missing	49.0	9.2	missing	missing	missing	missing
4	1	36.0	66.7	50	1	missing	0.952857	0.96443	missing	0	missing	32.0	8.5	missing	missing	missing	missing
5	1	48.0	66.7	50	1	missing	0.952857	0.96443	missing	0	missing	26.0	6.4	missing	missing	missing	missing

Step 3: Read the merged dataframe into a Pumas Population object.

# Convert resulting DataFrame to a Pumas Population object
pdpop_dataframe = read_pumas(
    pd_dataframe;
    id = :ID,
    time = :TIME,
    amt = :AMOUNT,
    cmt = :CMT,
    evid = :EVID,
    covariates = [:CL, :Vc, :tabs, :lags_Depot],
    observations = [:conc, :pca],
)

Population
  Subjects: 31
  Covariates: CL, Vc, tabs, lags_Depot
  Observations: conc, pca

Step 1: Extract the individual coefficients from the PK fit using icoef.

icoefs = icoef(pk_fit) # returns a vector of named tuples
dcs = dosecontrol(pk_fit) # returns a vector of named tuples

31-element Vector{Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}}:
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["1"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["2"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["3"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["4"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["5"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["6"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["7"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["8"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["9"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["11"], lags_Depot = [0.9261075051832584]), :left)
 ⋮
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["24"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["25"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["26"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["27"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["28"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["29"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["30"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["31"], lags_Depot = [0.9261075051832584]), :left)
 Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{id::Vector{String}, lags_Depot::Vector{Float64}}, Symbol}([0.0], (id = ["32"], lags_Depot = [0.9261075051832584]), :left)

Step 2: Merge the individual coefficients with the original population.

pdpop_subject_constructor_icoef = map(zip(icoefs, pkpd_pop)) do (row, pdsubj)
    Subject(pdsubj, covariates = row)
end

pdpop_subject_constructor_icoef_dcs =
    map(zip(dcs, pdpop_subject_constructor_icoef)) do (dc, pdsubj)
        Subject(pdsubj, covariates = dc)
    end

Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL, CL, Vc, tabs, Ka, lags_Depot
  Observations: conc, pca

First Step - Adding Individual PK Parameters:

The first step (pdpop_subject_constructor_icoef) adds the individual PK parameters (like clearance, volume, etc.) that were estimated from the PK model fit. These parameters come from icoef(pk_fit) which returns individual estimates for each subject.

Second Step - Adding Dose Control Parameters:

The second step (pdpop_subject_constructor_icoef_dcs) adds the dose control parameters (like lag times) that were also estimated from the PK model. These parameters come from dosecontrol(pk_fit) and are specifically related to how the drug dosing should be handled.

Step 1: Extract the individual coefficients from the PK fit (same as the DataFrames method).

Step 2: Extract the individual coefficients from the PK fit.

# Perform inspection on FittedPumasModel
warf_pk_insp = inspect(pk_fit);
# Convert to a DataFrame
df_inspect = DataFrame(warf_pk_insp);
# Obtain the individual PK parameters
icoef_dataframe = unique(df_inspect[!, [:id, :time, :CL, :Vc, :tabs, :lags_Depot]], :id)

first(icoef_dataframe, 5)

[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.

5×6 DataFrame

Row	id	time	CL	Vc	tabs	lags_Depot
	String	Float64	Float64?	Float64?	Float64?	Float64?
1	1	0.0	0.112039	7.67823	0.379683	0.926108
2	11	0.0	0.131106	7.7971	0.313889	0.926108
3	12	0.0	0.134624	6.83003	2.6345	0.926108
4	13	0.0	0.174219	8.78166	0.829443	0.926108
5	14	0.0	0.152349	6.15482	0.209459	0.926108

Step 3: Use the Subject constructor to create a new dataset for the PD model.

# We need to convert the dataframe to a vector of named tuples that can be mapped over and passed to the `Subject` constructor
icoef_nt = Tables.rowtable(icoef_dataframe)

pdpop_subject_constructor_icoef_dataframe = map(zip(icoef_nt, pkpd_pop)) do (row, pdsubj)
    Subject(
        pdsubj,
        covariates = (
            CL = row.CL,
            Vc = row.Vc,
            tabs = row.tabs,
            lags_Depot = row.lags_Depot,
        ),
    )
end

Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL, CL, Vc, tabs, lags_Depot
  Observations: conc, pca

Let’s break down this important code pattern that’s commonly used in sequential PK/PD modeling:

Converting Data Format

icoef_nt = Tables.rowtable(icoef_dataframe)

This converts our dataframe of individual PK parameters into named tuples
A named tuple is like a labeled container, for example:

(ID = "1", CL = 0.134, Vc = 8.11, tabs = 0.523, lags_Depot = 0.1)

Creating New Subjects

pdpop_subject_constructor_icoef_dataframe =
    map(zip(icoef_nt, pdpop_dataframe)) do (row, pdsubj)
        Subject(
            pdsubj,
            covariates = (
                CL = row.CL,
                Vc = row.Vc,
                tabs = row.tabs,
                lags_Depot = row.lags_Depot,
            ),
        )
    end

zip pairs each subject’s PK parameters with their original data
For each pair, we create a new subject that contains:
- All their original data (pdsubj)
- Their individual PK parameters as covariates

Real-World Analogy: Think of it like updating patient files:

Original file: Contains dosing and observation data
Update: Adding individual lab results (PK parameters) to each patient’s file

Before and After Example:

# Original Subject
Original_Subject = Subject(id = "1", observations = (conc = [10, 8, 6]), time = [0, 1, 2])

# After Adding PK Parameters
New_Subject = Subject(
    id = "1",
    observations = (conc = [10, 8, 6]),
    time = [0, 1, 2],
    covariates = (iCL = 0.134, iV = 8.11, itabs = 0.523, ilag = 0.1),
)

This pattern allows us to carry forward our individual PK parameter estimates into the PD modeling step while maintaining all original data.

3.2 Re-Defining the PD Model with the New Dataset


seq_warfarin_pkpd_model = @model begin

    @metadata begin
        desc = "Sequential Warfarin PK/PD model"
        timeu = u"hr"
    end

    @param begin
        # PD parameters
        """
        Baseline
        """
        pop_e0 ∈ RealDomain(lower = 0.0, init = 100.0)
        """
        Emax
        """
        pop_emax ∈ RealDomain(init = -1.0)
        """
        EC50
        """
        pop_c50 ∈ RealDomain(lower = 0.0, init = 1.0)
        """
        Turnover
        """
        pop_tover ∈ RealDomain(lower = 0.0, init = 14.0)
        # Inter-individual variability
        """
          - Ωe0
          - Ωemax
          - Ωec50
          - Ωturn
        """
        pd_Ω ∈ PDiagDomain([0.01, 0.01, 0.01, 0.01])
        # Residual variability
        """
        Proportional residual error for drug concentration
        """
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        """
        Additive residual error for drug concentration (mg/L)
        """
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
        """
        Additive error for PCA
        """
        σ_fx ∈ RealDomain(lower = 0.0, init = 0.01)
    end

    @random begin
        # mean = 0, covariance = pd_Ω
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates begin
        """
        Individual Clearance
        """
        CL
        """
        Individual Volume
        """
        Vc
        """
        Individual Absorption time
        """
        tabs
        """
        Individual Lag time
        """
        lags_Depot
    end

    @pre begin
        # PK
        Ka = log(2) / tabs
        # PD
        e0 = pop_e0 * exp(pd_η[1])
        emax = pop_emax * exp(pd_η[2])
        c50 = pop_c50 * exp(pd_η[3])
        tover = pop_tover * exp(pd_η[4])
        kout = log(2) / tover
        rin = e0 * kout
    end

    @dosecontrol begin
        lags = (Depot = lags_Depot,)
    end

    @init begin
        Turnover = e0
    end

    @vars begin
        cp := Central / Vc
        ratein := Ka * Depot
        pd := 1 + emax * cp / (c50 + cp)
    end

    @dynamics begin
        Depot' = -ratein
        Central' = ratein - CL * cp
        Turnover' = rin * pd - kout * Turnover
    end

    @derived begin
        """
        Warfarin Concentration (mg/L)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
        """
        PCA
        """
        pca ~ @. Normal(Turnover, σ_fx)
    end
end

PumasModel
  Parameters: pop_e0, pop_emax, pop_c50, pop_tover, pd_Ω, σ_prop, σ_add, σ_fx
  Random effects: pd_η
  Covariates: CL, Vc, tabs, lags_Depot
  Dynamical system variables: Depot, Central, Turnover
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca
  Observed: conc, pca

3.3 Fitting the PD Model (Sequential)

The loglikehood for the PD model given the two datasets derived above can be compared to ensure that both ways of deriving the PKPD data are equivalent.

ll_dataframe = loglikelihood(
    seq_warfarin_pkpd_model,
    pdpop_dataframe,
    init_params(seq_warfarin_pkpd_model),
    FOCE(),
)

-4.794218571769937e6

ll_subject_constructor_icoef_dcs = loglikelihood(
    seq_warfarin_pkpd_model,
    pdpop_subject_constructor_icoef_dcs,
    init_params(seq_warfarin_pkpd_model),
    FOCE(),
)

-4.794218571769935e6

Now we can check that the loglikelihoods are the same using both methods of deriving the PKPD data.

ll_dataframe ≈ ll_subject_constructor_icoef_dcs

true

Below is the PD model fit using the Subject constructor.

pdpop_subject_constructor_fit = fit(
    seq_warfarin_pkpd_model,
    pdpop_subject_constructor_icoef_dcs,
    init_params(seq_warfarin_pkpd_model),
    FOCE(),
);

FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:               Nonlinear ODE
Solver(s):(OrdinaryDiffEq.Vern7,OrdinaryDiffEq.Rodas5P)

Log-likelihood value:                   -894.33446
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0             11
Observation records:         Active        Missing
    conc:                       239             47
    pca:                        225             61
    Total:                      464            108

--------------------------
              Estimate
--------------------------
pop_e0        96.771
pop_emax      -1.0553
pop_c50        1.4644
pop_tover     13.852
pd_Ω₁,₁        0.0029729
pd_Ω₂,₂        9.6894e-10
pd_Ω₃,₃        0.17797
pd_Ω₄,₄        0.015551
σ_prop         0.0786
σ_add          0.28503
σ_fx           3.6777
--------------------------

We will not explore the results further here, but the PD model fit using the Subject constructor is complete.

In summary, the sequential approach facilitates a simpler conceptual model, but the direct communication of uncertainty from the PK parameters into the PD model is lost.

4 Simultaneous PK/PD Modeling

For a simultaneous approach, a joint model that handles both PK and PD data at once is built. Conceptually, the concentration variable is shared between PK and PD within the same model code, so all parameters are estimated in a single fit.

4.1 Defining the Joint PK/PD Model

The joint PK/PD model for sequential fitting has already been defined above in warfarin_pkpd_model.

PumasModel
  Parameters: pop_CL, pop_V, pop_tabs, pop_lag, pop_e0, pop_emax, pop_c50, pop_tover, pk_Ω, pd_Ω, σ_prop, σ_add, σ_fx
  Random effects: pk_η, pd_η
  Covariates: FSZV, FSZCL
  Dynamical system variables: Depot, Central, Turnover
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca
  Observed: conc, pca

There are two observed variables in this model:

conc for the observed Warfarin concentrations,
pca for the observed PD endpoints.

4.2 Building a Joint Dataset

A single dataset that contains both PK and PD data, possibly at different times, is needed. Pumas can handle that by allowing the definition of multiple observations while constructing a Population:

pkpd_pop = read_pumas(
    warfarin_data_wide;
    id = :ID,
    time = :TIME,
    amt = :AMOUNT,
    cmt = :CMT,
    evid = :EVID,
    covariates = [:SEX, :WEIGHT, :FSZV, :FSZCL],
    observations = [:conc, :pca],
)

Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL
  Observations: conc, pca

Here, conc holds the PK observations at certain times, of which may be missing at the PD observation times.
pca is present or may be missing at the PK observation times and has real values at the PD times.

4.3 Fitting the PK/PD Model (Simultaneous)

Given the joint dataset, the model can be fit using the same syntax as before:

simultaneous_pkpd_fit =
    fit(warfarin_pkpd_model, pkpd_pop, init_params(warfarin_pkpd_model), Pumas.FOCE());

FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:               Nonlinear ODE
Solver(s):(OrdinaryDiffEq.Vern7,OrdinaryDiffEq.Rodas5P)

Log-likelihood value:                   -1011.5354
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0             18
Observation records:         Active        Missing
    conc:                       239             47
    pca:                        225             61
    Total:                      464            108

-------------------------
              Estimate
-------------------------
pop_CL         0.13137
pop_V          8.0255
pop_tabs       0.42635
pop_lag        0.92336
pop_e0        96.688
pop_emax      -1.0628
pop_c50        1.5016
pop_tover     14.055
pk_Ω₁,₁        0.050744
pk_Ω₂,₂        0.021666
pk_Ω₃,₃        0.97884
pd_Ω₁,₁        0.0028828
pd_Ω₂,₂        7.3726e-9
pd_Ω₃,₃        0.15184
pd_Ω₄,₄        0.011913
σ_prop         0.088292
σ_add          0.36362
σ_fx           3.5336
-------------------------

That’s it! A single fit that estimates both PK and PD parameters, taking into account the full joint likelihood, is now available.

5 Comparison of Sequential and Simultaneous Approaches

We can compare the parameter estimates from the sequential and simultaneous approaches.

compare_estimates(;
    Sequential = pdpop_subject_constructor_fit,
    Simultaneous = simultaneous_pkpd_fit,
)

18×3 DataFrame

Row	parameter	Sequential	Simultaneous
	String	Float64?	Float64?
1	pop_e0	96.7709	96.6879
2	pop_emax	-1.05527	-1.06281
3	pop_c50	1.46439	1.50157
4	pop_tover	13.8517	14.0553
5	pd_Ω₁,₁	0.0029729	0.00288284
6	pd_Ω₂,₂	9.68939e-10	7.37255e-9
7	pd_Ω₃,₃	0.177971	0.151835
8	pd_Ω₄,₄	0.0155508	0.011913
9	σ_prop	0.0785999	0.088292
10	σ_add	0.285034	0.363621
11	σ_fx	3.67768	3.53356
12	pop_CL	missing	0.131369
13	pop_V	missing	8.02554
14	pop_tabs	missing	0.426354
15	pop_lag	missing	0.923356
16	pk_Ω₁,₁	missing	0.0507439
17	pk_Ω₂,₂	missing	0.0216655
18	pk_Ω₃,₃	missing	0.978843

6 Advantages and Disadvantages (Revisited)

By now, the key difference is clear:

Sequential:
1. Fit PK model.
2. Generate predicted concentrations.
3. Input predicted concentrations (as data) into the PD model.
Advantage: The PD model is simpler to interpret and fits quickly because it treats the Warfarin concentration as a known input. However, the direct communication of uncertainty from the PK parameters into the PD model is lost.
Simultaneous:
1. Fit PK and PD data in one step.
2. This approach provides the correct propagation of parameter uncertainty and potentially more robust estimates, but can be more complex to set up and can be slower or more sensitive to initial guesses.

Preparing the PKPD Data

Sequential approach: Both the DataFrame and the Subject constructor methods are easy to use, with the later being more flexible.
Simultaneous approach: The joint dataset can be prepared in a more flexible way, allowing for more complex data structures.

7 Conclusion

In this tutorial, the following was covered:

Defining a Warfarin PK model in Pumas.
Fitting that PK model to concentration data.
Performing sequential PK/PD modeling by using predicted PK concentrations as an input in a separate PD model.
Performing a simultaneous PK/PD approach by coding a joint model that handles PK and PD data in a single fit.

Reuse

CC BY-SA 4.0