using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using Random
using DataFrames
using DataFramesMeta
using Statistics
using StatsBase
TMDD Part 5: Practical Modeling Workflows
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:
- Part 1: Understanding TMDD - Conceptual foundations and parameter meaning 1b. Part 1b: TMDD Dynamics Through the Four Phases - Deep dive into phase characterization
- Part 2: The Full TMDD Model - Complete mechanistic model implementation
- Part 3: TMDD Approximations - QE, QSS, and Michaelis-Menten models
- Part 4: Special Cases and Extensions - Advanced TMDD scenarios
- 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
4 Workflow Overview
A typical TMDD modeling workflow includes:
- Data Exploration: Identify TMDD signatures
- Model Selection Strategy: Choose starting model
- Parameter Initialization: Set reasonable starting values
- Model Fitting: Estimate parameters
- Diagnostics: Evaluate model performance
- Model Comparison: Select best model
- 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 ~ @. Normal(CONC, 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 | η_1 | η_2 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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.00729686 | 2.16417 | 0.04 | 0.008 | 1.0 | 0.002 | 0.195684 | -0.326576 | 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",
),
)
Interpreting 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"),
)
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
endApparent 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
TMDD 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:
- TMDD is likely (dose-dependent half-life, non-superimposable profiles)
- Only ligand data available (no receptor measurements)
- Multiple doses (can estimate saturable kinetics)
Recommended approach:
- Start with Michaelis-Menten (MM) model - most parsimonious
- If MM is insufficient, try QSS model
- 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 ~ @. Normal(CONC, abs(CONC) * σ_prop)
end
endPumasModel
Parameters: tvCl, tvVc, tvVmax, tvKM, Ω, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Ltot
Dynamical system type: Nonlinear ODE
Derived: DV
Observed: DV8.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 ~ @. Normal(CONC, abs(CONC) * σ_prop)
end
endPumasModel
Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKSS, tvKint, Ω, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Ltot, Rtot
Dynamical system type: Nonlinear ODE
Derived: DV
Observed: DV9 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 nM9.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 nM10 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)4010.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.368318e+06 6.818270e+07 * time: 0.03497004508972168 1 3.518952e+03 5.651342e+03 * time: 1.953261137008667 2 3.518492e+03 5.650059e+03 * time: 2.023620128631592 3 1.688161e+03 1.778167e+03 * time: 2.17631196975708 4 1.281965e+03 9.247181e+02 * time: 2.285907030105591 5 1.059735e+03 5.487672e+02 * time: 2.3907880783081055 6 9.918427e+02 3.708877e+02 * time: 2.4756550788879395 7 9.734374e+02 2.824128e+02 * time: 2.5559160709381104 8 9.697935e+02 2.483885e+02 * time: 2.6356489658355713 9 9.682474e+02 2.399039e+02 * time: 2.724885940551758 10 9.680247e+02 2.399039e+02 * time: 3.043221950531006
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: NoXChange
Log-likelihood value: -968.02471
-------------------
Estimate
-------------------
tvCl 0.0060765
tvVc 4.9148
tvVmax 0.0035624
tvKM 2.4162
Ω₁,₁ 0.099895
Ω₂,₂ 0.099895
σ_prop 0.55258
-------------------# View parameter estimates
coef(fit_mm)(tvCl = 0.006076500237624084,
tvVc = 4.914778140154867,
tvVmax = 0.003562360343387881,
tvKM = 2.416216600577859,
Ω = [0.09989514614744037 0.0; 0.0 0.09989514614744024],
σ_prop = 0.5525786775900324,)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.112329e+02 5.936495e+02 * time: 1.8835067749023438e-5 1 7.723399e+02 3.883102e+02 * time: 0.7276849746704102 2 4.785962e+02 5.893334e+02 * time: 0.91444993019104 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/SciMLBase/JKXkh/src/integrator_interface.jl:589 3 3.239864e+02 3.910262e+02 * time: 2.393272876739502 4 2.025913e+02 1.580006e+03 * time: 2.5071158409118652 5 1.574160e+02 1.499932e+03 * time: 2.679142951965332 6 7.691523e+01 1.474377e+03 * time: 2.974167823791504 7 3.715658e+01 1.050838e+03 * time: 3.189605951309204 8 5.524321e+00 8.070437e+02 * time: 3.3860578536987305 9 -1.196141e+01 1.005780e+02 * time: 3.5283379554748535 10 -1.376596e+01 1.628545e+02 * time: 3.629326820373535 11 -2.234130e+01 1.181731e+02 * time: 3.7306790351867676 12 -2.674890e+01 3.329203e+01 * time: 3.8299028873443604 13 -4.190372e+01 8.838240e+01 * time: 3.945840835571289 14 -4.493484e+01 6.897034e+01 * time: 4.17343282699585 15 -4.803752e+01 5.016509e+01 * time: 4.380928993225098 16 -5.221539e+01 3.128154e+01 * time: 4.570791006088257 17 -5.271125e+01 5.705201e+01 * time: 4.680341958999634 18 -5.313534e+01 3.111548e+01 * time: 4.787869930267334 19 -5.316540e+01 2.102529e+01 * time: 4.907060861587524 20 -5.325000e+01 3.645377e+00 * time: 5.0111730098724365 21 -5.326419e+01 1.790435e+00 * time: 5.111699819564819 22 -5.332542e+01 1.041046e+01 * time: 5.215303897857666 23 -5.345118e+01 2.230396e+01 * time: 5.332389831542969 24 -5.368080e+01 3.319482e+01 * time: 5.439038038253784 25 -5.398219e+01 3.242593e+01 * time: 5.542776823043823 26 -5.432283e+01 1.362342e+01 * time: 5.657037973403931 27 -5.445277e+01 3.608022e+00 * time: 5.762430906295776 28 -5.445792e+01 8.339472e-01 * time: 5.85956883430481 29 -5.445807e+01 7.150681e-01 * time: 5.958624839782715 30 -5.445807e+01 1.116407e-01 * time: 6.062453985214233 31 -5.445808e+01 2.508574e-02 * time: 6.150542974472046 32 -5.445808e+01 1.502113e-02 * time: 6.235167980194092 33 -5.445808e+01 1.151705e-03 * time: 6.314935922622681 34 -5.445809e+01 1.109092e-04 * time: 6.393751859664917
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: GradientNorm
Log-likelihood value: 54.458091
-------------------
Estimate
-------------------
tvCl 0.0043988
tvVc 2.9515
tvKsyn 0.048684
tvKdeg 0.0094121
tvKSS 1.0279
tvKint 0.0021149
Ω₁,₁ 0.06194
Ω₂,₂ 0.06194
σ_prop 0.09821
-------------------# View parameter estimates
coef(fit_qss)(tvCl = 0.0043987699952257805,
tvVc = 2.951502555364329,
tvKsyn = 0.048683575682659436,
tvKdeg = 0.009412105350249863,
tvKSS = 1.0278877629344312,
tvKint = 0.002114910250396442,
Ω = [0.06194034248867061 0.0; 0.0 0.06194034248795285],
σ_prop = 0.09821010203217118,)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
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.
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,))
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")
endModel Comparison:
==================================================
Metric MM Model QSS Model
--------------------------------------------------
Log-likelihood -968.0 54.5
AIC 1950.0 -90.9
BIC 1979.8 -52.6
Parameters 6 8
==================================================
✓ QSS model preferred based on AIC12.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 |
Model 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
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
14 Workflow Summary
The complete TMDD modeling workflow:
Exploratory Analysis
- Plot dose-normalized concentrations
- Check for non-superimposable profiles
- Estimate apparent half-life by dose
Model Strategy
- Start with simplest model (MM)
- Progress to QSS if mechanistic insights needed
- Use full model only with receptor data
Parameter Initialization
- Use exploratory analysis to inform initial values
- Leverage simpler model fits for complex models
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
Model Selection
- Compare AIC/BIC
- Consider scientific objectives
- Balance parsimony with mechanism
Simulation
- Predict therapeutic scenarios
- Evaluate target engagement
- Support dose selection decisions
Key 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
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