Building a Compartmental PK Model

1 Introduction

This tutorial walks through the process of building a pharmacokinetic (PK) model of the combined PKPD model for Warfarin that was introduced in the previous tutorial. It is assumed that the two-compartment model with first-order absorption and absorption lag, inter-individual variability on clearance (\(CL)\), central volume of distribution (\(V_c\)), and absorption half-life (\(t_{abs}\)), allometric scaling on PK parameters, and a combined proportional and additive residual error model has already been established. Conceptually, the structural model can be summarized as follows:

1. Depot (absorption compartment)
2. Central (systemic)
3. Peripheral (distribution)

2 Breakdown of the Model

2.1 Model Structure

A schematic representation might look like this:

• Depot is where the drug is administered (e.g., via oral dosing).
• Central is the primary systemic compartment with clearance \(CL\).
• Peripheral is the secondary distribution compartment, exchanging drug with the central compartment via the inter-compartmental clearance \(Q\).
• \(Ka\) is the first-order absorption rate constant from the Depot to the Central compartment.

Allometric (size-based) scaling is often applied to parameters such as clearance (\(CL\)) and volume (\(V\)) to account for patient size differences or other covariate information. In this tutorial, scaling factors provided by covariates (e.g., FSZCL for clearance and FSZV for volume) will serve that role.

2.2 Writing the Differential Equations

Given the two-compartment PK structure and first-order absorption, the ordinary differential equations (ODEs) can be written as follows:

\[ \begin{aligned} \frac{d(\text{Depot})}{dt} &= - K_a \times \text{Depot}, \\[6pt] \frac{d(\text{Central})}{dt} &= K_a \times \text{Depot} \;-\; \frac{CL}{V_c} \times \text{Central} \;-\; \frac{Q}{V_c} \times \text{Central} \;+\; \frac{Q}{V_p} \times \text{Peripheral}, \\[6pt] \frac{d(\text{Peripheral})}{dt} &= \frac{Q}{V_c} \times \text{Central} \;-\; \frac{Q}{V_p} \times \text{Peripheral}. \end{aligned} \]

Here:

• \(K_a = \log(2)/t_{abs}\) is derived from an absorption half-life \(t_{abs}\).
• \(CL\) (clearance), \(V_c\) (central volume of distribution), \(Q\) (inter-compartmental clearance), and \(V_p\) (peripheral volume of distribution) can be made functions of covariates and random effects.
• Allometric scaling is handled by multiplying each base population parameter by the appropriate factor (e.g., FSZCL or FSZV) and exponentiating random effects if desired.

2.3 Converting to a @dynamics Block

In Pumas, the above system of ODEs can be written directly in the @dynamics block. For example:

@dynamics begin
    Depot' = -Ka * Depot
    Central' = Ka * Depot - (CL / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
    Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
end

This block automatically interprets each line as one of the state equations (e.g. \(\frac{d(\text{state})}{dt} = \dots\)).

2.4 Using the Coefficients in the @pre Block

Before these equations are solved, the model needs to specify how each parameter is computed from the population values, covariates, and random effects. This is typically performed in a @pre block. Each parameter in the ODE system (e.g., \(CL\), \(V_c\), \(Q\), \(V_p\), \(K_a\)) is expressed in terms of:

• Population parameters (e.g., pop_CL, pop_Vc, pop_Q, pop_Vp, pop_tabs)
• Covariates (e.g., FSZCL, FSZV)
• Random effects (e.g., exp(pk_η[i]))

An example set of algebraic expressions is shown below:

@pre begin
    # Example of allometric scaling using covariates:
    CL = pop_CL * FSZCL * exp(pk_η[1])
    Vc = pop_Vc * FSZV * exp(pk_η[2])

    # Inter-compartmental parameters without random effects in this tutorial
    Q = pop_Q * FSZCL
    Vp = pop_Vp * FSZV

    # Absorption half-life and rate constant
    tabs = pop_tabs * exp(pk_η[3])
    Ka = log(2) / tabs
end

In this illustration:

• FSZCL and FSZV are the covariate-based scaling factors.
• pop_CL, pop_Vc, pop_Q, and pop_Vp are the typical population values.
• pk_η[1], pk_η[2], and pk_η[3] are the individual random effects for \(CL\), \(V_c\), and \(t_{abs}\), respectively (in this example, \(Q\) and \(V_p\) have no inter-individual variability).

2.5 Specifying the Lag Time in the @dosecontrol Block

The lag time is specified in the @dosecontrol block. For example:

@dosecontrol begin
    lags = (Depot = pop_lag * exp(lag_η),)
end

2.6 Specifying Population Parameters in the @param Block

These population parameters must be declared so that the modeling environment recognizes them, along with any initial guesses, bounds, or descriptions. Below is a @param block that focuses on the PK parameters:

@param begin
    """
    Clearance (L/h/70kg)
    """
    pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)

    """
    Central Volume (L/70kg)
    """
    pop_Vc ∈ RealDomain(lower = 0.0, init = 8.11)

    """
    Inter-compartmental Clearance (L/h/70kg)
    """
    pop_Q ∈ RealDomain(lower = 0.0, init = 0.5)

    """
    Peripheral Volume (L/70kg)
    """
    pop_Vp ∈ RealDomain(lower = 0.0, init = 20.0)

    """
    Absorption time (h)
    """
    pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)

    """
    Lag time (h)
    """
    pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)

    # Inter-individual variability (diagonal covariance for CL, Vc, tabs)
    """
      - ΩCL
      - ΩVc
      - ΩTabs
    """
    pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])

    """
    Ωlag
    """
    lag_ω ∈ RealDomain(lower = 0.0, init = 0.1)

    # Residual variability for concentration
    """
    Proportional error
    """
    σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)

    """
    Additive error
    """
    σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
end

• Each parameter has a lower bound (e.g., 0.0) and an initial estimate (init) used during model fitting.
• pk_Ω captures the variances for the three random effects on \(CL\), \(V_c\), and \(t_{abs}\).
• lag_ω captures the standard deviation for the random effect on the lag time.

2.7 Specifying the Random Effects in the @random Block

The random effects are drawn from specified distributions in the @random block. For instance:

@random begin
    # mean = 0, variances = pk_Ω
    pk_η ~ MvNormal(pk_Ω)

    # mean = 0, standard deviation = lag_ω
    lag_η ~ Normal(0.0, lag_ω)
end

• pk_η is a vector of random effects for the PK parameters (CL, Vc, tabs, etc.).
• lag_η is a single random effect for the lag time.

During model simulation or estimation, each individual in the data set will have its own sample from these distributions, thereby introducing inter-individual variability.

2.8 Constructing the @derived Block for Predicted Concentration

The final step is to specify the model’s predictions for observed data—in this case, the drug concentration in plasma. The predicted concentration (cp) is calculated by dividing the mass in the central compartment by the central volume of distribution. An associated @vars block might look like this:

@vars begin
    cp := Central / Vc
end

This tells the model framework how to calculate cp from the states and parameters.

A residual error model is then applied to link cp to the observed concentration. For example:

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

• cp is the predicted concentration in the central compartment.
• σ_prop and σ_add capture proportional and additive measurement error, respectively.

2.9 Putting It All Together

Below is a concise example model that incorporates the steps above. All PD components from the original warfarin model have been removed, leaving only the two-compartment PK system (with a depot for absorption), IIV on \(CL\), \(V_c\), and \(t_{abs}\), allometric scaling factors (FSZCL, FSZV), and a combined proportional and additive residual variability for concentrations:

warfarin_pk_model = @model begin

    @metadata begin
        desc = "Warfarin 2-compartment PK model (PD removed)"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        pop_Vc ∈ RealDomain(lower = 0.0, init = 8.11)
        pop_Q ∈ RealDomain(lower = 0.0, init = 0.5)
        pop_Vp ∈ RealDomain(lower = 0.0, init = 20.0)
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
        pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)

        # Inter-individual variability
        pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])
        lag_ω ∈ RealDomain(lower = 0.0, init = 0.1)

        # Residual variability
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
    end

    @random begin
        # mean = 0, covariance = pk_Ω
        pk_η ~ MvNormal(pk_Ω)
        # mean = 0, standard deviation = lag_ω
        lag_η ~ Normal(0.0, lag_ω)
    end

    @covariates begin
        # Scaling factors
        FSZCL
        FSZV
    end

    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])
        Q = FSZCL * pop_Q
        Vp = FSZV * pop_Vp

        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
    end

    @dosecontrol begin
        lags = (Depot = pop_lag * exp(lag_η),)
    end

    @vars begin
        cp := Central / Vc
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' =
            Ka * Depot - (CL / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
        Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
    end

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

1. The @param block lists all population-level parameters.
2. The @random block defines how random effects are distributed.
3. The @covariates block provides covariate inputs for allometric scaling.
4. The @pre block transforms population parameters and random effects into individual-level parameters (\(CL, V_c, Q, V_p, K_a\)).
5. The @dosecontrol block defines the lag time for the absorption compartment.
6. The @vars block defines the predicted concentration.
7. The @dynamics block contains the ODEs for the two-compartment PK model.
8. The @derived block generates the predicted concentration conc with a specified residual error model.

Through these steps, the tutorial demonstrates how to move from known compartmental equations to a fully specified model that can be used for fitting or simulation in a population PK context.