PD Model Selection and Comparison

1 Learning Objectives

By the end of this tutorial, you will be able to:

• Apply a systematic approach to PD model selection
• Use exploratory data analysis to identify model type
• Generate and interpret hysteresis plots
• Compare competing models statistically
• Avoid common pitfalls in PD modeling
• Make scientifically justified model choices

2 Introduction

Selecting the appropriate PD model is crucial for accurate characterization of drug-effect relationships (Gabrielsson and Weiner 2016; Danhof et al. 2007). This tutorial provides a systematic framework for choosing between the four major PD model paradigms:

1. Direct Response: Immediate equilibrium
2. Effect Compartment: Delayed equilibrium (hysteresis)
3. Indirect Response: Turnover-mediated effects
4. K-PD: Irreversible binding mechanisms

flowchart TD
    Start[PK-PD Data] --> EDA[Exploratory Data Analysis]

    EDA --> Q1{Plot Effect vs.<br/>Concentration}
    Q1 --> |Single curve| Direct[Consider<br/>Direct Response]
    Q1 --> |Loop observed| Q2{Hysteresis<br/>direction?}

    Q2 --> |Counter-clockwise| EC[Consider<br/>Effect Compartment]
    Q2 --> |Clockwise| Q3{Tolerance or<br/>indirect mechanism?}

    Q3 --> |Tolerance| Tol[Tolerance Model]
    Q3 --> |Indirect| IDR[Consider IDR]

    EDA --> Q4{Effect persists<br/>after drug clears?}
    Q4 --> |Moderately| IDR
    Q4 --> |Dramatically<br/>days/weeks| KPD[Consider K-PD]

    EDA --> Q5{Known mechanism<br/>of action?}
    Q5 --> |Receptor binding| Direct
    Q5 --> |Enzyme turnover| IDR
    Q5 --> |Covalent binding| KPD

    style Direct fill:#e1f5fe
    style EC fill:#fff3e0
    style IDR fill:#e8f5e9
    style KPD fill:#fce4ec

Figure 1: Decision framework for PD model selection

3 Step 1: Exploratory Data Analysis

Before fitting any model, thorough exploratory analysis is essential.

3.1 Plotting Effect vs. Time

# Generate example data with different characteristics
# We'll create datasets that represent each model type

# Common PK component
pk_model = @model begin
    @param begin
        pop_Ka ∈ RealDomain(; lower = 0)
        pop_CL ∈ RealDomain(; lower = 0)
        pop_Vc ∈ RealDomain(; lower = 0)
    end
    @pre begin
        Ka = pop_Ka
        CL = pop_CL
        Vc = pop_Vc
    end
    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end
    @vars begin
        Cp = Central / Vc
    end
end

pk_params = (pop_Ka = 1.0, pop_CL = 0.5, pop_Vc = 20.0)
dosing = DosageRegimen(100, time = 0, cmt = 1)
subject = Subject(id = 1, events = dosing)
sim_pk = simobs(pk_model, subject, pk_params, obstimes = 0:0.5:24, simulate_error = false)

SimulatedObservations
  Simulated variables: Cp
  Time: 0.0:0.5:24.0

3.2 The Hysteresis Plot

The hysteresis plot (effect vs. concentration) is the most informative diagnostic:

# Convert simulation to DataFrame
pk_df = DataFrame(sim_pk)

# Generate three scenarios
E0, Emax, EC50 = 20.0, 80.0, 2.5
conc = pk_df.Cp
times = pk_df.time

# Direct response: single curve
effect_direct = E0 .+ (Emax .* conc) ./ (EC50 .+ conc)

# Counter-clockwise hysteresis (effect compartment)
# Simulate delayed concentration
ce = zeros(length(conc))
ke0 = 0.3
dt = 0.5
for i = 2:length(ce)
    ce[i] = ce[i-1] + dt * ke0 * (conc[i] - ce[i-1])
end
effect_ccw = E0 .+ (Emax .* ce) ./ (EC50 .+ ce)

# Clockwise hysteresis (tolerance)
# Effect decreases over time at same concentration
tolerance = 1.0 .- 0.5 .* (1 .- exp.(-0.1 .* times))
effect_cw = (E0 .+ (Emax .* conc) ./ (EC50 .+ conc)) .* tolerance

# Create DataFrames for each scenario
direct_df = DataFrame(conc = conc, effect = effect_direct, type = "direct")
ccw_df = DataFrame(conc = conc, effect = effect_ccw, type = "ccw")
cw_df = DataFrame(conc = conc, effect = effect_cw, type = "cw")

# Plot
fig = Figure(size = (900, 400))

# Direct response
direct_plt =
    data(direct_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :steelblue) +
        visual(Scatter, markersize = 8, color = :steelblue)
    )
draw!(fig[1, 1], direct_plt; axis = (title = "No Hysteresis\n(Direct Response)",))

# Counter-clockwise
ccw_plt =
    data(ccw_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :coral) +
        visual(Scatter, markersize = 8, color = :coral)
    )
ax2 = Axis(
    fig[1, 2],
    xlabel = "Concentration",
    ylabel = "Effect",
    title = "Counter-Clockwise\n(Effect Compartment)",
)
draw!(ax2, ccw_plt)
arrows!(
    ax2,
    [conc[8]],
    [effect_ccw[8]],
    [conc[12] - conc[8]],
    [effect_ccw[12] - effect_ccw[8]],
    arrowsize = 15,
    color = :red,
)

# Clockwise
cw_plt =
    data(cw_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :purple) +
        visual(Scatter, markersize = 8, color = :purple)
    )
ax3 = Axis(
    fig[1, 3],
    xlabel = "Concentration",
    ylabel = "Effect",
    title = "Clockwise\n(Tolerance/IDR)",
)
draw!(ax3, cw_plt)
arrows!(
    ax3,
    [conc[8]],
    [effect_cw[8]],
    [conc[4] - conc[8]],
    [effect_cw[4] - effect_cw[8]],
    arrowsize = 15,
    color = :red,
)

fig

Figure 2: Different hysteresis patterns and their model implications

3.3 Key Diagnostic Patterns

Observation	Interpretation	Candidate Model
No loop (single curve)	Immediate equilibrium	Direct Response
Counter-clockwise loop	Effect lags concentration	Effect Compartment
Clockwise loop	Tolerance or indirect mechanism	IDR or Tolerance
Effect >> drug duration	Irreversible or slow turnover	K-PD or IDR

4 Step 2: Mechanism-Based Selection

When pharmacological mechanism is known, use it to guide model selection (Mould and Upton 2012; Upton and Mould 2014):

Receptor Binding

Drug binds reversibly to receptor

• Rapid equilibration → Direct Response
• Slow equilibration → Effect Compartment

Examples: Beta-blockers, opioids, benzodiazepines

Enzyme Inhibition

Drug inhibits enzyme reversibly

• Fast on/off → Direct Response
• Product turnover involved → Indirect Response

Examples: Statins (HMG-CoA reductase), ACE inhibitors

Signal Transduction

Drug triggers downstream cascade

• Synthesis/degradation involved → Indirect Response
• Second messenger turnover → Indirect Response

Examples: Corticosteroids, growth factors

Covalent Binding

Drug binds irreversibly to target

• Target turnover governs recovery → K-PD

Examples: Aspirin, PPIs, covalent kinase inhibitors

5 Step 3: Statistical Model Comparison

5.1 Information Criteria

For comparing non-nested models, use information criteria:

AIC = -2 \cdot \log(L) + 2k

BIC = -2 \cdot \log(L) + k \cdot \log(n)

Where:

• L = maximum likelihood
• k = number of parameters
• n = number of observations

Interpretation

Lower AIC/BIC is better. AIC favors model fit while BIC penalizes complexity more heavily. A difference < 2 means models are essentially equivalent; 2-10 indicates moderate evidence for the better model; > 10 indicates strong evidence.

5.2 Likelihood Ratio Test

For nested models (one is a special case of another):

LRT = -2 \cdot (\log L_{reduced} - \log L_{full})

The LRT follows a chi-squared distribution with degrees of freedom equal to the difference in parameter count.

6 Case Study: Warfarin Anticoagulation

Let’s work through a real example comparing models for warfarin’s anticoagulant effect.

# Load warfarin data
warfarin_data = dataset("pumas/warfarin_pumas")

# Process data
warfarin_df = @chain warfarin_data begin
    @rtransform begin
        :id = string(:id)
        :FSZV = :wtbl / 70
        :FSZCL = (:wtbl / 70)^0.75
    end
end

# Create population for PKPD
pkpd_pop = read_pumas(
    warfarin_df;
    id = :id,
    time = :time,
    amt = :amt,
    cmt = :cmt,
    evid = :evid,
    covariates = [:sex, :wtbl, :FSZV, :FSZCL],
    observations = [:conc, :pca],
)

Population
  Subjects: 32
  Covariates: sex, wtbl, FSZV, FSZCL
  Observations: conc, pca

6.1 Model 1: Direct Response

warf_direct = @model begin
    @metadata begin
        desc = "Warfarin Direct Response"
        timeu = u"hr"
    end

    @param begin
        pop_CL ∈ RealDomain(; lower = 0, init = 0.13)
        pop_Vc ∈ RealDomain(; lower = 0, init = 8.0)
        pop_tabs ∈ RealDomain(; lower = 0, init = 0.5)
        pop_lag ∈ RealDomain(; lower = 0, init = 0.1)
        pop_E0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_Emax ∈ RealDomain(; init = -50.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 1.0)
        pk_Ω ∈ PDiagDomain([0.04, 0.04, 0.04])
        pd_Ω ∈ PDiagDomain([0.04])
        σ_prop ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add ∈ RealDomain(; lower = 0, init = 0.1)
        σ_pca ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates FSZCL FSZV

    @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
        E0 = pop_E0 * exp(pd_η[1])
        Emax = pop_Emax
        EC50 = pop_EC50
    end

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

    @dynamics Depots1Central1

    @vars begin
        cp = Central / Vc
        pca_pred = E0 + Emax * cp / (EC50 + cp)
    end

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

PumasModel
  Parameters: pop_CL, pop_Vc, pop_tabs, pop_lag, pop_E0, pop_Emax, pop_EC50, pk_Ω, pd_Ω, σ_prop, σ_add, σ_pca
  Random effects: pk_η, pd_η
  Covariates: FSZCL, FSZV
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: conc, pca, cp, pca_pred
  Observed: conc, pca, cp, pca_pred

6.2 Model 2: Indirect Response (Type III)

warf_idr = @model begin
    @metadata begin
        desc = "Warfarin IDR Type III"
        timeu = u"hr"
    end

    @param begin
        pop_CL ∈ RealDomain(; lower = 0, init = 0.13)
        pop_Vc ∈ RealDomain(; lower = 0, init = 8.0)
        pop_tabs ∈ RealDomain(; lower = 0, init = 0.5)
        pop_lag ∈ RealDomain(; lower = 0, init = 0.1)
        pop_E0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_Emax ∈ RealDomain(; init = 1.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 1.0)
        pop_turnover ∈ RealDomain(; lower = 0, init = 14.0)
        pk_Ω ∈ PDiagDomain([0.04, 0.04, 0.04])
        pd_Ω ∈ PDiagDomain([0.04])
        σ_prop ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add ∈ RealDomain(; lower = 0, init = 0.1)
        σ_pca ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates FSZCL FSZV

    @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
        E0 = pop_E0 * exp(pd_η[1])
        Emax = pop_Emax
        EC50 = pop_EC50
        kout = log(2) / pop_turnover
        kin0 = E0 * kout
    end

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

    @init begin
        PCA = E0
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - CL / Vc * Central
        PCA' = kin - kout * PCA
    end

    @vars begin
        # Plasma concentration
        cp = Central / Vc
        # Drug stimulation of production
        stim = Emax * cp / (EC50 + cp)
        kin = kin0 * (1 + stim)
        # PCA for plotting
        pca_pred = PCA
    end

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

PumasModel
  Parameters: pop_CL, pop_Vc, pop_tabs, pop_lag, pop_E0, pop_Emax, pop_EC50, pop_turnover, pk_Ω, pd_Ω, σ_prop, σ_add, σ_pca
  Random effects: pk_η, pd_η
  Covariates: FSZCL, FSZV
  Dynamical system variables: Depot, Central, PCA
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca, cp, stim, kin, pca_pred
  Observed: conc, pca, cp, stim, kin, pca_pred

6.3 Fit Both Models

# Fit direct response
fit_direct = fit(warf_direct, pkpd_pop, init_params(warf_direct), FOCE())

# Fit IDR
fit_idr = fit(warf_idr, pkpd_pop, init_params(warf_idr), FOCE())

6.4 Compare Models

# AIC comparison
aic_direct = aic(fit_direct)
aic_idr = aic(fit_idr)

println("AIC Comparison:")
println("  Direct Response: $(round(aic_direct, digits=1))")
println("  IDR Model:       $(round(aic_idr, digits=1))")
println("  Difference:      $(round(aic_direct - aic_idr, digits=1))")

AIC Comparison:
  Direct Response: 2437.4
  IDR Model:       2293.8
  Difference:      143.6

# BIC comparison
bic_direct = bic(fit_direct)
bic_idr = bic(fit_idr)

println("\nBIC Comparison:")
println("  Direct Response: $(round(bic_direct, digits=1))")
println("  IDR Model:       $(round(bic_idr, digits=1))")
println("  Difference:      $(round(bic_direct - bic_idr, digits=1))")

BIC Comparison:
  Direct Response: 2495.9
  IDR Model:       2356.5
  Difference:      139.5

6.5 Visual Comparison

# Inspect both models
insp_direct = inspect(fit_direct)
insp_idr = inspect(fit_idr)

[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.
[ 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

# GOF for direct model
goodness_of_fit(insp_direct; observations = [:pca])

Figure 3: Goodness-of-fit comparison for PCA observations

goodness_of_fit(insp_idr; observations = [:pca])

Figure 4: Goodness-of-fit for IDR model

7 Common Pitfalls

Pitfall 1: Overparameterization

Adding parameters without data support leads to:

• Poor parameter precision (wide confidence intervals)
• Model instability
• Overfitting

Solution: Start simple, add complexity only when data supports it.

Pitfall 2: Ignoring Mechanism

Choosing models based purely on fit statistics without considering:

• Known pharmacology
• Biological plausibility
• Clinical interpretation

Solution: Use mechanism to inform model selection, then validate with data.

Pitfall 3: Sparse PD Data

With limited PD observations:

• Hysteresis may not be detectable
• Turnover parameters poorly estimated
• Complex models not identifiable

Solution: Match model complexity to data richness.

8 Decision Checklist

Use this checklist when selecting a PD model:

1. Exploratory Analysis

   • Plot effect vs. time for all subjects
   • Create hysteresis plots (effect vs. concentration)
   • Assess variability in response

2. Mechanistic Considerations

   • Review drug’s mechanism of action
   • Identify target type (receptor, enzyme, etc.)
   • Consider known pharmacology

3. Model Fitting

   • Start with simplest appropriate model
   • Assess parameter identifiability
   • Check convergence and stability

4. Model Comparison

   • Compare nested models with LRT
   • Compare non-nested models with AIC/BIC
   • Evaluate goodness-of-fit diagnostics

5. Validation

   • Visual predictive checks
   • Residual analysis
   • Clinical plausibility of parameters

9 Summary

In this tutorial, we covered:

1. Systematic approach: EDA → Mechanism → Statistical comparison
2. Hysteresis analysis: Key diagnostic for model selection
3. Information criteria: AIC/BIC for model comparison
4. Case study: Warfarin model comparison
5. Common pitfalls: Overparameterization, ignoring mechanism, sparse data

9.1 Key Takeaways

• Always start with exploratory analysis
• Let mechanism inform model choice
• Use statistics to compare, not choose
• Simpler models are often better
• Parameter interpretability matters

9.2 Quick Reference Table

Data Pattern	Mechanism	Recommended Model
No hysteresis	Rapid receptor binding	Direct Response
Counter-clockwise loop	Distribution delay	Effect Compartment
Effect outlasts concentration	Turnover process	Indirect Response
Effect >> drug half-life	Irreversible binding	K-PD
Clockwise loop	Tolerance/adaptation	Tolerance Model

References

Danhof, M., J. de Jongh, E. C. De Lange, O. Della Pasqua, B. A. Ploeger, and R. A. Voskuyl. 2007. “Mechanism-Based Pharmacokinetic-Pharmacodynamic Modeling: Biophase Distribution, Receptor Theory, and Dynamical Systems Analysis.” Annual Review of Pharmacology and Toxicology 47: 357–400.

Gabrielsson, J., and D. Weiner. 2016. Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications. 5th ed. Swedish Pharmaceutical Press.

Mould, D. R., and R. N. Upton. 2012. “Basic Concepts in Population Modeling, Simulation, and Model-Based Drug Development.” CPT: Pharmacometrics & Systems Pharmacology 1 (9): e6.

Upton, R. N., and D. R. Mould. 2014. “Basic Concepts in Population Modeling, Simulation, and Model-Based Drug Development: Part 3—Introduction to Pharmacodynamic Modeling Methods.” CPT: Pharmacometrics & Systems Pharmacology 3 (1): e88.