TMDD Part 5: Practical Modeling Workflows

Author

Vijay Ivaturi

1 Introduction

This is Part 5 of our comprehensive TMDD tutorial series. Here we bring together everything we’ve learned into practical end-to-end workflows for TMDD model development, from exploratory analysis through model fitting and selection.

Tutorial Series Overview:

  1. Part 1: Understanding TMDD - Conceptual foundations and parameter meaning 1b. Part 1b: TMDD Dynamics Through the Four Phases - Deep dive into phase characterization
  2. Part 2: The Full TMDD Model - Complete mechanistic model implementation
  3. Part 3: TMDD Approximations - QE, QSS, and Michaelis-Menten models
  4. Part 4: Special Cases and Extensions - Advanced TMDD scenarios
  5. Part 5: Practical Workflows (this tutorial) - Fitting, diagnostics, and model selection

2 Learning Objectives

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

  • Perform exploratory data analysis to identify TMDD signatures
  • Select appropriate initial parameter values for TMDD models
  • Fit TMDD models using a stepwise approach (MM to QSS to Full)
  • Generate diagnostic plots specific to TMDD
  • Compare models using information criteria and scientific reasoning
  • Simulate therapeutic scenarios for decision-making

3 Libraries

using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using Random
using DataFrames
using DataFramesMeta
using Statistics
using StatsBase

4 Workflow Overview

A typical TMDD modeling workflow includes:

  1. Data Exploration: Identify TMDD signatures
  2. Model Selection Strategy: Choose starting model
  3. Parameter Initialization: Set reasonable starting values
  4. Model Fitting: Estimate parameters
  5. Diagnostics: Evaluate model performance
  6. Model Comparison: Select best model
  7. Simulation and Prediction: Use model for decision-making

5 Generating Simulated Data

For this tutorial, we’ll generate simulated data from a known TMDD model to demonstrate the workflow. In practice, you would use your own clinical or preclinical data.

Show/Hide Code 
# Define the "true" TMDD model for data generation
true_model = @model begin
    @param begin
        tvCl ∈ RealDomain(lower = 0)
        tvVc ∈ RealDomain(lower = 0)
        tvKsyn ∈ RealDomain(lower = 0)
        tvKdeg ∈ RealDomain(lower = 0)
        tvKSS ∈ RealDomain(lower = 0)
        tvKint ∈ RealDomain(lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_prop ∈ RealDomain(lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Cl = tvCl * exp(η[1])
        Vc = tvVc * exp(η[2])
        Ksyn = tvKsyn
        Kdeg = tvKdeg
        KSS = tvKSS
        Kint = tvKint
    end

    @dosecontrol begin
        bioav = (Ltot = 1 / tvVc,)
    end

    @init begin
        Rtot = Ksyn / Kdeg
    end

    @vars begin
        Kel := Cl / Vc
        R0 := Ksyn / Kdeg
        L := 0.5 * ((Ltot - Rtot - KSS) + sqrt((Ltot - Rtot - KSS)^2 + 4 * KSS * Ltot))
        LR := Ltot - L
        R := Rtot - LR
    end

    @dynamics begin
        Ltot' = -Kel * L - Kint * LR
        Rtot' = Ksyn - Kdeg * R - Kint * LR
    end

    @derived begin
        CONC := @. L
        DV ~ @. ProportionalNormal(CONC, σ_prop)
    end
end

# True parameter values based on typical mAb parameters
true_params = (
    tvCl = 0.006,      # ~0.15 L/day
    tvVc = 3.0,
    tvKsyn = 0.04,     # ~1 nmol/day
    tvKdeg = 0.008,    # ~0.2 day⁻¹
    tvKSS = 1.0,       # 1 nM
    tvKint = 0.002,    # ~0.04 day⁻¹
    Ω = Diagonal([0.04, 0.04]),  # 20% IIV on Cl and Vc
    σ_prop = 0.1,  # 10% proportional error
)

# Create population: 40 subjects across 4 dose groups
Random.seed!(123)

doses = [5.0, 15.0, 50.0, 150.0]  # nmol (wide range for clear TMDD visualization)
n_per_group = 10

# Observation times (typical sparse sampling for mAbs)
obs_times = [0.5, 1, 2, 4, 8, 24, 48, 72, 168, 336, 504, 672, 1008]

sim_population = Subject[]
for (dose_idx, dose) in enumerate(doses)
    for subj = 1:n_per_group
        id = (dose_idx - 1) * n_per_group + subj
        push!(
            sim_population,
            Subject(
                id = id,
                events = DosageRegimen(dose, time = 0, cmt = 1),
                covariates = (DOSE = dose, WT = 70.0),
                observations = (DV = nothing,),
            ),
        )
    end
end

# Simulate data
sim_data = simobs(true_model, sim_population, true_params, obstimes = obs_times)

# Convert to DataFrame
df = DataFrame(sim_data)
@rtransform! df :LLOQ = 0.01;
@rsubset! df :DV > :LLOQ;  # Remove BLQ observations
@rtransform! df :route = "IV";
first(df, 10)
10×27 DataFrame
Row id time DV evid bioav_Ltot amt cmt rate duration ss ii route DOSE WT tad dosenum Ltot Rtot η₁ η₂ Cl Vc Ksyn Kdeg KSS Kint LLOQ
String Float64 Float64? Int64? Float64? Float64? Symbol? Float64? Float64? Int8? Float64? String Float64? Float64? Float64? Int64? Float64? Float64? Float64? Float64? Float64? Float64? Float64? Float64? Float64? Float64? Float64
1 1 0.5 0.369569 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 0.5 1 1.66476 5.00392 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
2 1 1.0 0.337308 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 1.0 1 1.66285 5.00783 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
3 1 2.0 0.346385 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 2.0 1 1.65904 5.01558 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
4 1 4.0 0.316171 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 4.0 1 1.65146 5.03086 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
5 1 8.0 0.307751 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 8.0 1 1.63641 5.06054 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
6 1 24.0 0.308507 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 24.0 1 1.5777 5.16811 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
7 1 48.0 0.2777 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 48.0 1 1.49391 5.29961 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
8 1 72.0 0.286412 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 72.0 1 1.41487 5.40087 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
9 1 168.0 0.187265 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 168.0 1 1.13998 5.59876 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01
10 1 336.0 0.153013 0 missing 0.0 missing 0.0 0.0 0 0.0 IV 5.0 70.0 336.0 1 0.78303 5.57511 0.195684 -0.326576 0.00729686 2.16417 0.04 0.008 1.0 0.002 0.01

6 Step 1: Exploratory Data Analysis

6.1 Dose-Normalized Concentration Plots

The hallmark of TMDD is non-superimposable dose-normalized profiles. Let’s check for this signature:

Show/Hide Code 
# Add dose information and calculate dose-normalized concentration
df_plot = @chain df begin
    @rtransform :DOSE_NORM = :DV / :DOSE
    @rsubset :time > 0
    @rtransform :Dose_Label = string(:DOSE) * " nmol"
end

plt =
    data(df_plot) *
    mapping(
        :time => "Time (hr)",
        :DOSE_NORM => "Dose-Normalized Concentration";
        color = :Dose_Label => "Dose",
    ) *
    visual(Scatter; markersize = 8, alpha = 0.5)

draw(
    plt;
    axis = (
        yscale = Makie.pseudolog10,
        title = "Exploratory Analysis: Dose-Normalized Profiles",
    ),
)

Dose-normalized concentration profiles. Non-superimposition at terminal phase suggests TMDD.
NoteInterpreting Dose-Normalized Plots

If profiles are superimposable, linear PK is likely. If lower doses show faster terminal decline (as we see here), TMDD is a likely explanation.

6.2 Mean Concentration by Dose

Show/Hide Code 
df_summary = @chain df begin
    @rsubset :time > 0
    groupby([:time, :DOSE])
    @combine begin
        :mean_DV = mean(:DV)
        :lower = quantile(:DV, 0.05)
        :upper = quantile(:DV, 0.95)
    end
    @rtransform :Dose_Label = string(:DOSE) * " nmol"
end

# Create line and band layers
plt_line =
    data(df_summary) *
    mapping(
        :time => "Time (hr)",
        :mean_DV => "Mean Concentration (nM)";
        color = :Dose_Label => "Dose",
    ) *
    visual(Lines; linewidth = 2)

plt_band =
    data(df_summary) *
    mapping(:time, :lower, :upper; color = :Dose_Label) *
    visual(Band; alpha = 0.2)

draw(
    plt_band + plt_line;
    axis = (yscale = Makie.pseudolog10, title = "Mean Concentration-Time Profiles"),
)

Mean concentration-time profiles by dose group with 90% confidence intervals.

6.3 Apparent Terminal Half-Life by Dose

Show/Hide Code 
# Calculate approximate apparent terminal half-life from terminal slope
# Using data from the terminal phase (last few time points)

println("Apparent Terminal Half-Life by Dose (estimated from terminal slope):")
for dose in doses
    # Get terminal phase data for this dose (last 3 time points with positive DV)
    df_dose = filter(row -> row.DOSE == dose && row.DV > 0, df)
    sort!(df_dose, :time)

    # Use data from later time points for terminal slope
    df_terminal = filter(row -> row.time >= 336, df_dose)  # After 2 weeks

    if nrow(df_terminal) >= 2
        # Calculate mean log-concentration at each time point
        df_mean = @chain df_terminal begin
            groupby(:time)
            @combine :mean_logDV = mean(log.(:DV))
        end

        if nrow(df_mean) >= 2
            # Linear regression on log scale to get terminal slope
            times = df_mean.time
            log_conc = df_mean.mean_logDV

            # Simple linear regression: slope = Σ(x-x̄)(y-ȳ) / Σ(x-x̄)²
            t_mean = mean(times)
            lc_mean = mean(log_conc)
            slope =
                sum((times .- t_mean) .* (log_conc .- lc_mean)) /
                sum((times .- t_mean) .^ 2)

            # Half-life = ln(2) / |slope|
            half_life = -log(2) / slope
            println("  $(dose) nmol: $(round(half_life, digits=1)) hr")
        end
    end
end
Apparent Terminal Half-Life by Dose (estimated from terminal slope):
  5.0 nmol: 341.5 hr
  15.0 nmol: 335.2 hr
  50.0 nmol: 292.2 hr
  150.0 nmol: 185.9 hr
ImportantTMDD Signature Confirmed

The increasing apparent half-life with dose is a classic TMDD signature. At low doses, target-mediated clearance dominates, leading to faster elimination.

7 Step 2: Model Selection Strategy

Based on our exploratory analysis:

  1. TMDD is likely (dose-dependent half-life, non-superimposable profiles)
  2. Only ligand data available (no receptor measurements)
  3. Multiple doses (can estimate saturable kinetics)

Recommended approach:

  1. Start with Michaelis-Menten (MM) model - most parsimonious
  2. If MM is insufficient, try QSS model
  3. Full model only if QSS fails and additional data become available

8 Step 3: Define Candidate Models

8.1 Michaelis-Menten Model

mm_model = @model begin
    @param begin
        tvCl ∈ RealDomain(lower = 0)
        tvVc ∈ RealDomain(lower = 0)
        tvVmax ∈ RealDomain(lower = 0)
        tvKM ∈ RealDomain(lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_prop ∈ RealDomain(lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Cl = tvCl * exp(η[1])
        Vc = tvVc * exp(η[2])
        Vmax = tvVmax
        KM = tvKM
    end

    @dosecontrol begin
        bioav = (Ltot = 1 / tvVc,)
    end

    @dynamics begin
        Ltot' = -(Cl / Vc) * Ltot - Vmax * Ltot / (KM + Ltot)
    end

    @derived begin
        CONC := @. Ltot
        DV ~ @. ProportionalNormal(CONC, σ_prop)
    end
end
PumasModel
  Parameters: tvCl, tvVc, tvVmax, tvKM, Ω, σ_prop
  Random effects: η
  Covariates:
  Dynamical system variables: Ltot
  Dynamical system type: Nonlinear ODE
  Derived: DV
  Observed: DV

8.2 Quasi-Steady-State (QSS) Model

qss_model = @model begin
    @param begin
        tvCl ∈ RealDomain(lower = 0)
        tvVc ∈ RealDomain(lower = 0)
        tvKsyn ∈ RealDomain(lower = 0)
        tvKdeg ∈ RealDomain(lower = 0)
        tvKSS ∈ RealDomain(lower = 0)
        tvKint ∈ RealDomain(lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_prop ∈ RealDomain(lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Cl = tvCl * exp(η[1])
        Vc = tvVc * exp(η[2])
        Ksyn = tvKsyn
        Kdeg = tvKdeg
        KSS = tvKSS
        Kint = tvKint
    end

    @dosecontrol begin
        bioav = (Ltot = 1 / tvVc,)
    end

    @init begin
        Rtot = Ksyn / Kdeg
    end

    @vars begin
        Kel := Cl / Vc
        R0 := Ksyn / Kdeg
        L := 0.5 * ((Ltot - Rtot - KSS) + sqrt((Ltot - Rtot - KSS)^2 + 4 * KSS * Ltot))
        LR := Ltot - L
        R := Rtot - LR
    end

    @dynamics begin
        Ltot' = -Kel * L - Kint * LR
        Rtot' = Ksyn - Kdeg * R - Kint * LR
    end

    @derived begin
        CONC := @. L
        DV ~ @. ProportionalNormal(CONC, σ_prop)
    end
end
PumasModel
  Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKSS, tvKint, Ω, σ_prop
  Random effects: η
  Covariates:
  Dynamical system variables: Ltot, Rtot
  Dynamical system type: Nonlinear ODE
  Derived: DV
  Observed: DV

9 Step 4: Parameter Initialization

Choosing good initial values is crucial for TMDD model convergence. Peletier & Gabrielsson (2012) provides guidance on extracting parameter estimates from different phases of the TMDD concentration-time profile.

9.1 Four-Phase Parameter Identification Strategy

Phase Observable Feature Initial Estimate Strategy
A Rapid initial decline \(k_{on}\) from duration (~\(1/(k_{on} \cdot C_0)\))
B Slope during saturation \(k_{el} = CL/V_c\) from linear portion
C Inflection point \(K_D\) or \(K_{SS}\) from transition concentration
D Terminal slope \(k_{int}\) from \(\lambda_z \approx k_{int}\)

Additionally:

  • Initial concentration (\(C_0\)) at high dose estimates \(V_c\)
  • Dose-dependence of terminal slope estimates \(R_0\) and binding parameters
  • Time to receptor recovery informs \(k_{syn}\) (since \(R_0 = k_{syn}/k_{deg}\))

9.2 Strategy for MM Model

# Initial estimates for MM model
# Based on exploratory analysis

# 1. Vc from early time points (assume rapid distribution)
#    Approximate from C0 = Dose/Vc
df_high_dose_early = filter(row -> row.DOSE == 150.0 && row.time < 1, df)
mean_c0_high_dose = mean(df_high_dose_early.DV)

Vc_init = 150.0 / mean_c0_high_dose  # High dose / C0

# 2. Cl from terminal slope at high dose (where TMDD is saturated)
#    Approximate: Cl ≈ Vc * kel (typical mAb ~0.15 L/day)
Cl_init = 0.006  # Start with typical mAb clearance

# 3. Vmax from difference in clearance at low vs high dose
#    Vmax = R0 * kint (for R0 ~ 5 nM, kint ~ 0.002 hr⁻¹)
Vmax_init = 0.01  # Estimate from dose-dependent clearance

# 4. KM from concentration where half-maximal saturation occurs
#    Typically near KSS for TMDD drugs
KM_init = 1.0

init_mm = (
    tvCl = Cl_init,
    tvVc = Vc_init,
    tvVmax = Vmax_init,
    tvKM = KM_init,
    Ω = Diagonal([0.1, 0.1]),
    σ_prop = 0.15,
)

println("MM Model Initial Estimates:")
println("  Cl: $(init_mm.tvCl) L/hr")
println("  Vc: $(round(init_mm.tvVc, digits=2)) L")
println("  Vmax: $(init_mm.tvVmax) nM/hr")
println("  KM: $(init_mm.tvKM) nM")
MM Model Initial Estimates:
  Cl: 0.006 L/hr
  Vc: 3.2 L
  Vmax: 0.01 nM/hr
  KM: 1.0 nM

9.3 Strategy for QSS Model

# Initial estimates for QSS model
# Based on typical mAb values from the literature

init_qss = (
    tvCl = 0.006,      # ~0.15 L/day
    tvVc = 3.0,
    tvKsyn = 0.04,     # ~1 nmol/day → R0 = 5 nM
    tvKdeg = 0.008,    # ~0.2 day⁻¹ (receptor t½ ~87 hr)
    tvKSS = 1.0,       # ≈ KM from MM model
    tvKint = 0.002,    # ~0.04 day⁻¹
    Ω = Diagonal([0.1, 0.1]),
    σ_prop = 0.15,
)

println("\nQSS Model Initial Estimates:")
println("  Cl: $(init_qss.tvCl) L/hr")
println("  Vc: $(init_qss.tvVc) L")
println("  Ksyn: $(init_qss.tvKsyn) nM/hr")
println("  Kdeg: $(init_qss.tvKdeg) hr⁻¹")
println("  KSS: $(init_qss.tvKSS) nM")
println("  Kint: $(init_qss.tvKint) hr⁻¹")
println("  Implied R0: $(init_qss.tvKsyn / init_qss.tvKdeg) nM")

QSS Model Initial Estimates:
  Cl: 0.006 L/hr
  Vc: 3.0 L
  Ksyn: 0.04 nM/hr
  Kdeg: 0.008 hr⁻¹
  KSS: 1.0 nM
  Kint: 0.002 hr⁻¹
  Implied R0: 5.0 nM

10 Step 5: Model Fitting

10.1 Prepare Data for Fitting

# Convert simulated observations to a fitting-ready population
# The simobs output contains the actual simulated observations
fitting_pop = Subject.(sim_data)

length(fitting_pop)
40

10.2 Fit Michaelis-Menten Model

# Fit MM model
fit_mm = fit(mm_model, fitting_pop, init_mm, 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.364434e+06     6.818388e+07
 * time: 0.0343778133392334
     1     3.519213e+03     5.651334e+03
 * time: 2.9081640243530273
     2     3.518594e+03     5.650052e+03
 * time: 2.959514856338501
     3     1.687761e+03     1.776548e+03
 * time: 3.013497829437256
     4     1.282388e+03     9.243676e+02
 * time: 3.0643489360809326
     5     1.061823e+03     5.486075e+02
 * time: 3.1156609058380127
     6     9.929534e+02     3.707387e+02
 * time: 3.16617488861084
     7     9.742209e+02     2.823710e+02
 * time: 3.21323299407959
     8     9.714914e+02     2.483856e+02
 * time: 3.2596359252929688
     9     9.713384e+02     2.397794e+02
 * time: 3.4109718799591064
    10     9.713352e+02     2.388137e+02
 * time: 3.4519150257110596
    11     9.713327e+02     2.383435e+02
 * time: 3.48842191696167
    12     9.713216e+02     2.369857e+02
 * time: 3.5252888202667236
    13     9.712970e+02     2.351086e+02
 * time: 3.561478853225708
    14     9.712283e+02     2.317302e+02
 * time: 3.5977249145507812
    15     9.710545e+02     2.261265e+02
 * time: 3.636586904525757
    16     9.706023e+02     2.162634e+02
 * time: 3.6779048442840576
    17     9.694623e+02     1.986852e+02
 * time: 3.71915602684021
    18     9.666983e+02     1.745125e+02
 * time: 3.8265528678894043
    19     9.607506e+02     1.721701e+02
 * time: 3.863579034805298
    20     9.512752e+02     1.464048e+02
 * time: 3.9047489166259766
    21     9.418000e+02     1.387710e+02
 * time: 3.9477698802948
    22     9.332831e+02     1.314877e+02
 * time: 3.990118980407715
    23     9.308423e+02     1.286175e+02
 * time: 4.032932996749878
    24     9.307221e+02     1.297247e+02
 * time: 4.072722911834717
    25     9.307202e+02     1.296801e+02
 * time: 4.113924980163574
    26     9.307197e+02     1.296728e+02
 * time: 4.179018020629883
    27     9.307155e+02     1.296276e+02
 * time: 4.214260816574097
    28     9.307074e+02     1.295651e+02
 * time: 4.252461910247803
    29     9.306832e+02     1.294229e+02
 * time: 4.288534879684448
    30     9.306227e+02     1.291367e+02
 * time: 4.330066919326782
    31     9.304610e+02     1.284860e+02
 * time: 4.370190858840942
    32     9.300401e+02     1.269775e+02
 * time: 4.410547971725464
    33     9.289415e+02     1.233412e+02
 * time: 4.478890895843506
    34     9.261505e+02     1.145870e+02
 * time: 4.526841878890991
    35     9.195523e+02     9.456213e+01
 * time: 4.5774829387664795
    36     9.071592e+02     1.349089e+02
 * time: 4.632100820541382
    37     8.943907e+02     1.195659e+02
 * time: 4.6943199634552
    38     8.886761e+02     5.422672e+01
 * time: 4.778235912322998
    39     8.872060e+02     2.710995e+01
 * time: 4.833794832229614
    40     8.871050e+02     2.471533e+01
 * time: 4.892765998840332
    41     8.871043e+02     2.449244e+01
 * time: 4.941339015960693
    42     8.871043e+02     2.448809e+01
 * time: 5.0060508251190186
    43     8.871041e+02     2.444573e+01
 * time: 5.050829887390137
    44     8.871037e+02     2.440328e+01
 * time: 5.098407983779907
    45     8.871026e+02     2.432028e+01
 * time: 5.146807909011841
    46     8.870999e+02     2.419252e+01
 * time: 5.192831993103027
    47     8.870925e+02     2.397320e+01
 * time: 5.263871908187866
    48     8.870735e+02     2.360382e+01
 * time: 5.313395023345947
    49     8.870238e+02     2.296024e+01
 * time: 5.363200902938843
    50     8.868957e+02     2.182429e+01
 * time: 5.417749881744385
    51     8.865706e+02     2.749364e+01
 * time: 5.471300840377808
    52     8.857736e+02     4.136772e+01
 * time: 5.543864011764526
    53     8.839389e+02     5.772643e+01
 * time: 5.5972490310668945
    54     8.802225e+02     6.921488e+01
 * time: 5.651752948760986
    55     8.750152e+02     5.837689e+01
 * time: 5.705978870391846
    56     8.711918e+02     2.182791e+01
 * time: 5.77829384803772
    57     8.697930e+02     1.671071e+01
 * time: 5.8291850090026855
    58     8.694603e+02     1.951235e+01
 * time: 5.876224994659424
    59     8.665494e+02     2.538943e+01
 * time: 5.926687002182007
    60     8.647183e+02     3.467844e+01
 * time: 6.058679819107056
    61     8.610116e+02     6.252763e+01
 * time: 6.189846992492676
    62     8.501078e+02     1.054226e+02
 * time: 6.233170032501221
    63     8.466983e+02     3.275473e+01
 * time: 6.331082820892334
    64     8.456780e+02     8.965742e+00
 * time: 6.365073919296265
    65     8.456140e+02     2.517938e+00
 * time: 6.4022979736328125
    66     8.456120e+02     3.061814e-01
 * time: 6.435933828353882
    67     8.456120e+02     2.937864e-02
 * time: 6.466768026351929
    68     8.456120e+02     4.546395e-03
 * time: 6.494282960891724
    69     8.456120e+02     4.546395e-03
 * time: 6.575270891189575
    70     8.456120e+02     7.674538e-03
 * time: 6.604362964630127
    71     8.456120e+02     6.998250e-03
 * time: 6.637022972106934
    72     8.456120e+02     6.313822e-03
 * time: 6.669155836105347
    73     8.456120e+02     5.682749e-03
 * time: 6.721380949020386
    74     8.456120e+02     5.677089e-03
 * time: 6.762728929519653
    75     8.456120e+02     5.671403e-03
 * time: 6.806579828262329
    76     8.456120e+02     5.665729e-03
 * time: 6.848994016647339
    77     8.456120e+02     5.660063e-03
 * time: 6.892773866653442
    78     8.456120e+02     5.654402e-03
 * time: 6.955157041549683
    79     8.456120e+02     5.648748e-03
 * time: 6.999406814575195
    80     8.456120e+02     5.648183e-03
 * time: 7.04712700843811
    81     8.456120e+02     5.647619e-03
 * time: 7.094743013381958
    82     8.456120e+02     5.647053e-03
 * time: 7.162590980529785
    83     8.456120e+02     5.646489e-03
 * time: 7.211717844009399
    84     8.456120e+02     5.645924e-03
 * time: 7.25959587097168
    85     8.456120e+02     5.645360e-03
 * time: 7.306959867477417
    86     8.456120e+02     5.644796e-03
 * time: 7.3766138553619385
    87     8.456120e+02     5.644739e-03
 * time: 7.431213855743408
    88     8.456120e+02     5.644683e-03
 * time: 7.484139919281006
    89     8.456120e+02     5.644678e-03
 * time: 7.563501834869385
    90     8.456120e+02     5.036251e-03
 * time: 7.595824956893921
    91     8.456120e+02     4.534151e-03
 * time: 7.629776954650879
    92     8.456120e+02     4.529631e-03
 * time: 7.672767877578735
    93     8.456120e+02     4.525089e-03
 * time: 7.715134859085083
    94     8.456120e+02     4.520562e-03
 * time: 7.777467966079712
    95     8.456120e+02     4.516041e-03
 * time: 7.821061849594116
    96     8.456120e+02     4.511525e-03
 * time: 7.864289999008179
    97     8.456120e+02     4.507013e-03
 * time: 7.90819787979126
    98     8.456120e+02     4.506562e-03
 * time: 7.976538896560669
    99     8.456120e+02     4.506112e-03
 * time: 8.024526834487915
   100     8.456120e+02     4.505661e-03
 * time: 8.073148012161255
   101     8.456120e+02     4.505209e-03
 * time: 8.120620965957642
   102     8.456120e+02     4.504760e-03
 * time: 8.1892409324646
   103     8.456120e+02     4.504309e-03
 * time: 8.237663984298706
   104     8.456120e+02     4.503859e-03
 * time: 8.286031007766724
   105     8.456120e+02     4.503408e-03
 * time: 8.333678960800171
   106     8.456120e+02     4.503363e-03
 * time: 8.408282041549683
   107     8.456120e+02     4.503318e-03
 * time: 8.462530851364136
   108     8.456120e+02     4.503274e-03
 * time: 8.515559911727905
   109     8.456120e+02     4.503228e-03
 * time: 8.59009599685669
   110     8.456120e+02     4.503183e-03
 * time: 8.64520788192749
   111     8.456120e+02     4.503138e-03
 * time: 8.699323892593384
   112     8.456120e+02     4.503134e-03
 * time: 8.757581949234009
   113     8.456120e+02     4.503129e-03
 * time: 8.839018821716309
   114     8.456120e+02     4.503124e-03
 * time: 8.899449825286865
   115     8.456120e+02     4.503121e-03
 * time: 8.957912921905518
   116     8.456120e+02     4.503115e-03
 * time: 9.037052869796753
   117     8.456120e+02     4.503111e-03
 * time: 9.095418930053711
   118     8.456120e+02     4.503106e-03
 * time: 9.153954029083252
   119     8.456120e+02     4.503102e-03
 * time: 9.21189284324646
   120     8.456120e+02     4.503097e-03
 * time: 9.290822982788086
   121     8.456120e+02     4.503096e-03
 * time: 9.355470895767212
   122     8.456120e+02     4.503096e-03
 * time: 9.42389988899231
   123     8.456120e+02     4.503097e-03
 * time: 9.524910926818848
   124     8.456120e+02     4.503097e-03
 * time: 9.612368822097778
   125     8.456120e+02     4.213909e-03
 * time: 9.65330696105957
   126     8.456120e+02     4.181613e-03
 * time: 9.713326930999756
   127     8.456120e+02     4.178367e-03
 * time: 9.755992889404297
   128     8.456120e+02     4.175133e-03
 * time: 9.801694869995117
   129     8.456120e+02     4.171904e-03
 * time: 9.845258951187134
   130     8.456120e+02     4.171580e-03
 * time: 9.917616844177246
   131     8.456120e+02     4.171258e-03
 * time: 9.967074871063232
   132     8.456120e+02     4.170934e-03
 * time: 10.018596887588501
   133     8.456120e+02     4.170611e-03
 * time: 10.068647861480713
   134     8.456120e+02     4.170287e-03
 * time: 10.141400814056396
   135     8.456120e+02     4.170255e-03
 * time: 10.19813084602356
   136     8.456120e+02     4.170221e-03
 * time: 10.253835916519165
   137     8.456120e+02     4.170218e-03
 * time: 10.315024852752686
   138     8.456120e+02     4.170215e-03
 * time: 10.398185968399048
   139     8.456120e+02     4.170212e-03
 * time: 10.460577011108398
   140     8.456120e+02     4.170211e-03
 * time: 10.548280000686646
   141     8.456120e+02     4.170212e-03
 * time: 10.623708963394165
   142     8.456120e+02     4.170212e-03
 * time: 10.712993860244751
   143     8.456120e+02     4.170212e-03
 * time: 10.800220966339111

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             40

Observation records:         Active        Missing
    DV:                         520              0
    Total:                      520              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:              NoObjectiveChange
Log-likelihood value:                   -845.61199

------------------------
         Estimate
------------------------
tvCl         0.018989
tvVc         5.5428
tvVmax       0.00024241
tvKM     47368.0
Ω₁,₁         0.0034195
Ω₂,₂         0.0034195
σ_prop       0.51214
------------------------
# View parameter estimates
coef(fit_mm)
(tvCl = 0.018989092348128882, tvVc = 5.542806138090067, tvVmax = 0.0002424139504615842, tvKM = 47367.99121761668, Ω = [0.003419537697620273 0.0; 0.0 0.003419537697683042], σ_prop = 0.512143657304441)

10.3 Fit QSS Model

# Use MM results to inform QSS initial values
mm_coefs = coef(fit_mm)

# Update QSS initial values based on MM fit
init_qss_updated = (
    tvCl = mm_coefs.tvCl,
    tvVc = mm_coefs.tvVc,
    tvKsyn = 0.04,      # ~1 nmol/day → R0 = 5 nM
    tvKdeg = 0.008,     # ~0.2 day⁻¹
    tvKSS = mm_coefs.tvKM,  # KM ≈ KSS
    tvKint = mm_coefs.tvVmax / 5,  # Vmax/R0 ≈ Kint (R0 ≈ 5)
    Ω = mm_coefs.Ω,
    σ_prop = mm_coefs.σ_prop,
)

# Fit QSS model
fit_qss = fit(qss_model, fitting_pop, init_qss_updated, 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     8.456120e+02     3.364426e-01
 * time: 1.6927719116210938e-5
     1     8.456120e+02     1.547311e-01
 * time: 1.3961570262908936
     2     8.456120e+02     1.218466e-01
 * time: 1.926259994506836
     3     8.456120e+02     1.230834e-01
 * time: 2.0652518272399902
     4     8.456120e+02     6.209838e-02
 * time: 3.639029026031494
     5     8.456120e+02     5.952754e-02
 * time: 3.731865882873535
     6     8.456120e+02     5.952754e-02
 * time: 3.936138868331909
     7     8.456120e+02     5.675769e-02
 * time: 4.040306806564331
     8     8.456120e+02     5.524358e-02
 * time: 4.129786014556885
     9     8.456120e+02     5.523216e-02
 * time: 4.22786283493042
    10     8.456120e+02     5.399013e-02
 * time: 4.378058910369873
    11     8.456120e+02     5.279127e-02
 * time: 4.4564208984375
    12     8.456120e+02     5.271078e-02
 * time: 4.5414488315582275
    13     8.456120e+02     5.264431e-02
 * time: 4.630644798278809
    14     8.456120e+02     5.257574e-02
 * time: 4.761714935302734
    15     8.456120e+02     5.195255e-02
 * time: 4.840817928314209
    16     8.456120e+02     5.137605e-02
 * time: 4.923983812332153
    17     8.456120e+02     5.137605e-02
 * time: 5.067312002182007

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             40

Observation records:         Active        Missing
    DV:                         520              0
    Total:                      520              0

Number of parameters:      Constant      Optimized
                                  0              9

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                      NoXChange
Log-likelihood value:                   -845.61198

-----------------------
         Estimate
-----------------------
tvCl         0.018989
tvVc         5.5422
tvKsyn       0.040002
tvKdeg       0.0079996
tvKSS    47364.0
tvKint       4.8483e-5
Ω₁,₁         0.0034197
Ω₂,₂         0.0034197
σ_prop       0.51215
-----------------------
# View parameter estimates
coef(fit_qss)
(tvCl = 0.018989246211922342, tvVc = 5.5421732034806634, tvKsyn = 0.04000154361767202, tvKdeg = 0.007999646915130998, tvKSS = 47364.222488029125, tvKint = 4.848263298119195e-5, Ω = [0.0034196679943480565 0.0; 0.0 0.0034196680760198717], σ_prop = 0.5121485169973713)

11 Step 6: Model Diagnostics

11.1 Goodness-of-Fit: Observations vs Predictions

Show/Hide Code 
# Generate predictions
pred_mm = predict(fit_mm)
pred_qss = predict(fit_qss)

# Extract predictions
df_mm_pred = DataFrame(pred_mm)
df_mm_pred.Model .= "MM Model"

df_qss_pred = DataFrame(pred_qss)
df_qss_pred.Model .= "QSS Model"

df_gof = vcat(
    select(df_mm_pred, :DV_pred, :DV, :Model),
    select(df_qss_pred, :DV_pred, :DV, :Model),
)

# Create scatter layer
plt =
    data(df_gof) *
    mapping(:DV_pred => "Population Predictions", :DV => "Observations"; layout = :Model) *
    visual(Scatter; markersize = 6, alpha = 0.4)

fg = draw(plt; axis = (aspect = 1,))

# Add line of identity to each panel
for i = 1:2
    ablines!(fg.figure[1, i], 0, 1; color = :black, linestyle = :dash)
end

fg

Observed vs population predictions for both models with line of identity.

11.2 Residual Diagnostics

Show/Hide Code 
# Calculate residuals
insp_mm = inspect(fit_mm)
insp_qss = inspect(fit_qss)

df_insp_mm = DataFrame(insp_mm)
df_insp_mm.Model .= "MM Model"

df_insp_qss = DataFrame(insp_qss)
df_insp_qss.Model .= "QSS Model"

df_resid = vcat(
    select(df_insp_mm, :time, :DV_wres, :Model),
    select(df_insp_qss, :time, :DV_wres, :Model),
)

# Create scatter layer
plt =
    data(df_resid) *
    mapping(:time => "Time (hr)", :DV_wres => "CWRES"; layout = :Model) *
    visual(Scatter; markersize = 6, alpha = 0.4)

fg = draw(plt)

# Add zero reference line to each panel
for i = 1:2
    hlines!(fg.figure[1, i], [0]; color = :gray, linestyle = :dash)
end

fg
[ 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.

Conditional weighted residuals (CWRES) vs time for both models.

11.3 Individual Fits

Show/Hide Code 
# Select one subject from each dose group
selected_ids = ["1", "11", "21", "31"]

# Create dose lookup from original data
dose_lookup = Dict(row.id => row.DOSE for row in eachrow(unique(select(df, :id, :DOSE))))

# Prepare data for individual fits
df_obs_select = filter(row -> row.id in selected_ids, df)
df_obs_select = @rtransform df_obs_select :Type = "Observed"
df_obs_select = @rtransform df_obs_select :Subject =
    "Subject " * string(:id) * " (" * string(:DOSE) * " nmol)"
df_obs_select = select(df_obs_select, :time, :DV, :Type, :Subject)
rename!(df_obs_select, :DV => :Conc)

# For inspect DataFrames, add DOSE from lookup
df_mm_select = filter(row -> row.id in selected_ids, df_insp_mm)
df_mm_select.DOSE = [dose_lookup[id] for id in df_mm_select.id]
df_mm_select = @rtransform df_mm_select :Type = "MM IPRED"
df_mm_select = @rtransform df_mm_select :Subject =
    "Subject " * string(:id) * " (" * string(:DOSE) * " nmol)"
df_mm_select = select(df_mm_select, :time, :DV_ipred, :Type, :Subject)
rename!(df_mm_select, :DV_ipred => :Conc)

df_qss_select = filter(row -> row.id in selected_ids, df_insp_qss)
df_qss_select.DOSE = [dose_lookup[id] for id in df_qss_select.id]
df_qss_select = @rtransform df_qss_select :Type = "QSS IPRED"
df_qss_select = @rtransform df_qss_select :Subject =
    "Subject " * string(:id) * " (" * string(:DOSE) * " nmol)"
df_qss_select = select(df_qss_select, :time, :DV_ipred, :Type, :Subject)
rename!(df_qss_select, :DV_ipred => :Conc)

df_indiv = vcat(df_obs_select, df_mm_select, df_qss_select)

# Create scatter for observed and lines for predictions
plt_obs =
    data(filter(row -> row.Type == "Observed", df_indiv)) *
    mapping(:time => "Time (hr)", :Conc => "Concentration (nM)"; layout = :Subject) *
    visual(Scatter; color = :black, markersize = 8)

plt_pred =
    data(filter(row -> row.Type != "Observed", df_indiv)) *
    mapping(
        :time => "Time (hr)",
        :Conc => "Concentration (nM)";
        color = :Type,
        layout = :Subject,
    ) *
    visual(Lines; linewidth = 2)

draw(plt_obs + plt_pred; axis = (yscale = Makie.pseudolog10,))

Individual fits for selected subjects across dose groups.

12 Step 7: Model Comparison

12.1 Information Criteria

# Compare models using AIC and BIC
println("Model Comparison:")
println("="^50)
println("Metric\t\tMM Model\tQSS Model")
println("-"^50)
println(
    "Log-likelihood\t$(round(loglikelihood(fit_mm), digits=1))\t\t$(round(loglikelihood(fit_qss), digits=1))",
)
println("AIC\t\t$(round(aic(fit_mm), digits=1))\t\t$(round(aic(fit_qss), digits=1))")
println("BIC\t\t$(round(bic(fit_mm), digits=1))\t\t$(round(bic(fit_qss), digits=1))")
println("Parameters\t$(length(coef(fit_mm)))\t\t$(length(coef(fit_qss)))")
println("="^50)

# Determine better model
if aic(fit_qss) < aic(fit_mm)
    println("\n✓ QSS model preferred based on AIC")
else
    println("\n✓ MM model preferred based on AIC")
end
Model Comparison:
==================================================
Metric      MM Model    QSS Model
--------------------------------------------------
Log-likelihood  -845.6      -845.6
AIC     1705.2      1709.2
BIC     1735.0      1747.5
Parameters  6       8
==================================================

✓ MM model preferred based on AIC

12.2 Scientific Considerations

Beyond statistical criteria, consider:

Aspect MM Model QSS Model
Receptor dynamics Cannot predict Can predict
Target engagement Cannot estimate Can estimate
Mechanistic interpretation Limited Moderate
Identifiability Good May need constraints
Extrapolation Risky Better mechanistic basis
TipModel Selection Recommendation

For this dataset:

  • If only PK predictions are needed, the MM model is adequate
  • If target engagement or receptor dynamics are needed, use the QSS model
  • The QSS model provides more mechanistic insights for dose selection

13 Step 8: Simulation for Decision-Making

13.1 Dose-Response Predictions

Using the fitted model, we can simulate therapeutic scenarios:

Show/Hide Code 
# Use the best-fit model for predictions
best_fit = fit_qss
best_model = qss_model
best_params = coef(best_fit)

# Simulation parameters (remove random effects for typical subject)
sim_params = (
    tvCl = best_params.tvCl,
    tvVc = best_params.tvVc,
    tvKsyn = best_params.tvKsyn,
    tvKdeg = best_params.tvKdeg,
    tvKSS = best_params.tvKSS,
    tvKint = best_params.tvKint,
    Ω = Diagonal([0.0, 0.0]),  # No IIV for typical subject
    σ_prop = 0.0,
)

# Define dosing scenarios (using nmol doses)
scenarios = [
    (dose = 15.0, interval = 168, label = "15 nmol Q1W"),
    (dose = 50.0, interval = 336, label = "50 nmol Q2W"),
    (dose = 150.0, interval = 672, label = "150 nmol Q4W"),
]

sim_times = 0:1:2016  # 12 weeks

# Simulate all scenarios and combine into DataFrame
df_scenarios = DataFrame()
for scenario in scenarios
    n_doses = ceil(Int, 2016 / scenario.interval) + 1
    regimen = DosageRegimen(
        scenario.dose,
        time = 0,
        cmt = 1,
        addl = n_doses - 1,
        ii = scenario.interval,
    )

    pop = [Subject(id = 1, events = regimen, observations = (DV = nothing,))]

    sim = simobs(qss_model, pop, sim_params, obstimes = sim_times)
    sim_df = DataFrame(sim)
    sim_df.Regimen .= scenario.label
    append!(df_scenarios, sim_df)
end

plt =
    data(df_scenarios) *
    mapping(
        :time => "Time (hr)",
        :DV => "Free Ligand Concentration (nM)";
        color = :Regimen,
    ) *
    visual(Lines; linewidth = 2)

fg = draw(
    plt;
    axis = (yscale = Makie.pseudolog10, title = "Simulated PK for Different Regimens"),
)

# Add KSS reference line
hlines!(fg.figure[1, 1], [best_params.tvKSS], color = :gray, linestyle = :dash)

fg

Simulated concentration-time profiles for proposed dosing regimens.

13.2 Target Suppression Predictions

Show/Hide Code 
# Compute receptor fraction from state variables after simulation
# Using the QSS relationship: R = Rtot - LR, where LR = Ltot - L

# Parameters for receptor calculation
R0 = sim_params.tvKsyn / sim_params.tvKdeg  # Baseline receptor
KSS = sim_params.tvKSS

# Calculate receptor fraction from df_scenarios
df_receptor = copy(df_scenarios)

# Compute free ligand L from QSS relationship
# Note: We use DV as proxy for free ligand concentration L
# For receptor, we need to work backwards from the state variables
# Since Ltot and Rtot aren't directly in output, we'll compute R_frac
# using the steady-state approximation at each time point

# At quasi-steady-state: L ≈ DV (free ligand)
# LR = R0 * L / (KSS + L) (from QSS assumption with constant R approximation)
# R = R0 - LR = R0 * KSS / (KSS + L)
# R_frac = R / R0 = KSS / (KSS + L)

df_receptor.L = df_receptor.DV
df_receptor.R_frac = KSS ./ (KSS .+ df_receptor.L)

plt =
    data(df_receptor) *
    mapping(
        :time => "Time (hr)",
        :R_frac => "Free Receptor (fraction of baseline)";
        color = :Regimen,
    ) *
    visual(Lines; linewidth = 2)

fg = draw(plt; axis = (title = "Predicted Target Suppression",))

# Add reference lines
hlines!(fg.figure[1, 1], [1.0], color = :gray, linestyle = :dash)
hlines!(fg.figure[1, 1], [0.1], color = :red, linestyle = :dot)

fg

Predicted receptor suppression (as fraction of baseline) for different regimens.

14 Workflow Summary

The complete TMDD modeling workflow:

  1. Exploratory Analysis

    • Plot dose-normalized concentrations
    • Check for non-superimposable profiles
    • Estimate apparent half-life by dose

  2. Model Strategy

    • Start with simplest model (MM)
    • Progress to QSS if mechanistic insights needed
    • Use full model only with receptor data

  3. Parameter Initialization

    • Use exploratory analysis to inform initial values
    • Leverage simpler model fits for complex models

  4. Fitting and Diagnostics

    • Check convergence
    • Evaluate residuals and GOF plots
    • Verify individual fits
    • Use four-phase analysis (Peletier & Gabrielsson, 2012) to verify model captures expected TMDD behavior

  5. Model Selection

    • Compare AIC/BIC
    • Consider scientific objectives
    • Balance parsimony with mechanism

  6. Simulation

    • Predict therapeutic scenarios
    • Evaluate target engagement
    • Support dose selection decisions
TipKey Takeaways
  • Always start with exploratory analysis to identify TMDD signatures
  • Begin with simple models and add complexity as needed
  • Good initial values are critical for convergence
  • Model selection should consider both statistics and scientific purpose
  • The final model should support the intended application

For comprehensive reviews of TMDD modeling approaches, see Gibiansky et al. (2008) and Dua et al. (2015). Population modeling concepts are discussed in Mould & Upton (2012).

15 References

Dua, P., Hawkins, E., & Graaf, P. H. van der. (2015). A tutorial on target-mediated drug disposition (TMDD) models. CPT: Pharmacometrics & Systems Pharmacology, 4(6), 324–337. https://doi.org/10.1002/psp4.41
Gibiansky, L., Gibiansky, E., Kakkar, T., & Ma, P. (2008). Approximations of the target-mediated drug disposition model and identifiability of model parameters. Journal of Pharmacokinetics and Pharmacodynamics, 35(5), 573–591. https://doi.org/10.1007/s10928-008-9102-8
Mould, D. R., & Upton, R. N. (2012). Basic concepts in population modeling, simulation, and model-based drug development. CPT: Pharmacometrics & Systems Pharmacology, 1, e6. https://doi.org/10.1038/psp.2012.4
Peletier, L. A., & Gabrielsson, J. (2012). Dynamics of target-mediated drug disposition: Characteristic profiles and parameter identification. Journal of Pharmacokinetics and Pharmacodynamics, 39(5), 429–451. https://doi.org/10.1007/s10928-012-9260-6

Reuse