using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using Random
using DataFrames
TMDD Part 3: Model Approximations (QE, QSS, Michaelis-Menten)
1 Introduction
This is Part 3 of our comprehensive TMDD tutorial series. Here we explore the practical approximations of the full TMDD model that are commonly used when the full model is either overparameterized or when data are insufficient for estimating all parameters.
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 (this tutorial) - QE, QSS, and Michaelis-Menten models
- Part 4: Special Cases and Extensions - Advanced TMDD scenarios
- Part 5: Practical Workflows - Fitting, diagnostics, and model selection
2 Learning Objectives
By the end of this tutorial, you will be able to:
- Explain why TMDD model simplifications are necessary
- Implement Quasi-Equilibrium (QE), Quasi-Steady-State (QSS), and Michaelis-Menten models in Pumas
- Understand the mathematical basis and assumptions behind each approximation
- Select the appropriate approximation based on your data and scientific context
- Compare model predictions across different approximation levels
3 Libraries
4 Why Simplify the Full TMDD Model?
The full TMDD model contains nine kinetic parameters: \(CL\), \(V_c\), \(Q\), \(V_p\), \(k_{syn}\), \(k_{deg}\), \(k_{on}\), \(k_{off}\), and \(k_{int}\). In practice, estimating all these parameters is often challenging because:
Parameter non-identifiability: With typical PK data (free ligand concentrations only), several parameters are structurally or practically non-identifiable.
Correlation between parameters: \(k_{on}\) and \(k_{off}\) are often highly correlated, making individual estimation unreliable.
Data limitations: Receptor and complex concentrations are rarely measured, limiting the information available for parameter estimation.
Numerical instability: The full model can have stiff ODEs that are computationally expensive and prone to convergence issues.
Gibiansky et al. (2008) provided critical insights into when simplifications are appropriate and how to select the right approximation.
5 The Approximation Hierarchy
The TMDD model approximations form a hierarchy from most mechanistic to most simplified:
Full TMDD → Quasi-Steady-State (QSS) → Quasi-Equilibrium (QE) → Michaelis-Menten (MM)
↑ ↑ ↑
Most general Assumes rapid Assumes receptor
binding kinetics concentration negligible
Each approximation makes additional assumptions that reduce model complexity but restrict applicability.
6 The Full TMDD Model (Reference)
For comparison, here is the one-compartment full TMDD model:
\[ \frac{dL}{dt} = -k_{el} \cdot L - k_{on} \cdot L \cdot R + k_{off} \cdot LR \]
\[ \frac{dR}{dt} = k_{syn} - k_{deg} \cdot R - k_{on} \cdot L \cdot R + k_{off} \cdot LR \]
\[ \frac{dLR}{dt} = k_{on} \cdot L \cdot R - k_{off} \cdot LR - k_{int} \cdot LR \]
Where:
- \(L\) = Free ligand concentration
- \(R\) = Free receptor concentration
- \(LR\) = Ligand-receptor complex concentration
- \(k_{el}\) = Linear elimination rate constant
- \(k_{on}\), \(k_{off}\) = Association and dissociation rate constants
- \(k_{syn}\), \(k_{deg}\) = Receptor synthesis and degradation rates
- \(k_{int}\) = Complex internalization rate
tmdd_full_model = @model begin
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKon ∈ RealDomain(lower = 0)
tvKoff ∈ RealDomain(lower = 0)
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
Ksyn = tvKsyn
Kdeg = tvKdeg
Kon = tvKon
Koff = tvKoff
Kint = tvKint
end
@dosecontrol begin
bioav = (Ligand = 1 / tvVc,)
end
@init begin
Receptor = Ksyn / Kdeg
end
@vars begin
Kel := Cl / Vc
R0 := Ksyn / Kdeg
Ltot := Ligand + Complex
Rtot := Receptor + Complex
end
@dynamics begin
Ligand' = -Kel * Ligand - Kon * Ligand * Receptor + Koff * Complex
Receptor' = Ksyn - Kdeg * Receptor - Kon * Ligand * Receptor + Koff * Complex
Complex' = Kon * Ligand * Receptor - Koff * Complex - Kint * Complex
end
@derived begin
Free_Ligand = @. Ligand
Free_Receptor = @. Receptor
Complex_Conc = @. Complex
Total_Ligand = @. Ltot
Total_Receptor = @. Rtot
end
endPumasModel
Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKon, tvKoff, tvKint
Random effects:
Covariates:
Dynamical system variables: Ligand, Receptor, Complex
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
Observed: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
7 Quasi-Equilibrium (QE) Approximation
7.1 Mathematical Basis
The QE approximation assumes that the binding and unbinding of ligand to receptor is much faster than all other processes (elimination, distribution, internalization). Mathematically, this means:
\[ k_{on} \cdot L \cdot R \approx k_{off} \cdot LR \]
This implies that the binding reaction is at equilibrium at all times, and we can define the equilibrium dissociation constant:
\[ K_D = \frac{k_{off}}{k_{on}} = \frac{L \cdot R}{LR} \]
7.2 When is QE Appropriate?
The QE approximation is valid when:
- \(k_{on}\) and \(k_{off}\) are both very large compared to \(k_{int}\), \(k_{deg}\), and \(k_{el}\)
- Specifically, \(k_{int} \ll k_{off}\) is a key requirement
QE assumes that the complex dissociates much faster than it is internalized. If \(k_{int}\) is comparable to or larger than \(k_{off}\), use QSS instead.
7.3 Deriving the QE Equations
From the equilibrium assumption and mass balance:
Total ligand: \(L_{tot} = L + LR\)
Total receptor: \(R_{tot} = R + LR\)
Equilibrium: \(K_D = \frac{L \cdot R}{LR}\)
Solving for free ligand \(L\) as a function of \(L_{tot}\) and \(R_{tot}\):
\[ L = \frac{1}{2}\left[(L_{tot} - R_{tot} - K_D) + \sqrt{(L_{tot} - R_{tot} - K_D)^2 + 4 \cdot K_D \cdot L_{tot}}\right] \]
The complex concentration is then:
\[ LR = L_{tot} - L \]
And free receptor:
\[ R = R_{tot} - LR \]
7.4 QE Model in Pumas
tmdd_qe_model = @model begin
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKD ∈ RealDomain(lower = 0) # KD replaces Kon/Koff
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
Ksyn = tvKsyn
Kdeg = tvKdeg
KD = tvKD
Kint = tvKint
end
@dosecontrol begin
bioav = (Ltot = 1 / tvVc,)
end
@init begin
Rtot = Ksyn / Kdeg
end
@vars begin
Kel := Cl / Vc
R0 := Ksyn / Kdeg
# Solve for free ligand from equilibrium
L := 0.5 * ((Ltot - Rtot - KD) + sqrt((Ltot - Rtot - KD)^2 + 4 * KD * Ltot))
# Complex and free receptor from mass balance
LR := Ltot - L
R := Rtot - LR
end
@dynamics begin
Ltot' = -Kel * L - Kint * LR
Rtot' = Ksyn - Kdeg * R - Kint * LR
end
@derived begin
Free_Ligand = @. L
Free_Receptor = @. R
Complex_Conc = @. LR
Total_Ligand = @. Ltot
Total_Receptor = @. Rtot
end
endPumasModel
Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKD, tvKint
Random effects:
Covariates:
Dynamical system variables: Ltot, Rtot
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
Observed: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
Notice that QE reduces the number of differential equations from 3 to 2, and replaces \(k_{on}\) and \(k_{off}\) with a single parameter \(K_D\). This improves identifiability when only ligand concentrations are measured.
8 Quasi-Steady-State (QSS) Approximation
8.1 Mathematical Basis
The QSS approximation, introduced by Gibiansky et al. (2008), is more general than QE. It assumes that the rate of complex formation equals the rate of complex removal:
\[ k_{on} \cdot L \cdot R = (k_{off} + k_{int}) \cdot LR \]
This leads to the steady-state constant:
\[ K_{SS} = \frac{k_{off} + k_{int}}{k_{on}} \]
Note that \(K_{SS} \geq K_D\) always, with equality only when \(k_{int} = 0\).
8.2 When is QSS Appropriate?
The QSS approximation is valid when:
- Binding kinetics are fast relative to other processes
- \(k_{int}\) may be significant (unlike QE)
- Both ligand and receptor concentrations are measured, OR
- There is a strong scientific rationale for the target biology
The key difference is in the steady-state constant:
- QE: \(K_D = k_{off} / k_{on}\) (assumes \(k_{int} \ll k_{off}\))
- QSS: \(K_{SS} = (k_{off} + k_{int}) / k_{on}\) (accounts for internalization)
QSS is more general and should be preferred when \(k_{int}\) is not negligible.
8.3 QSS Model in Pumas
tmdd_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) # KSS replaces Kon/Koff/Kint combination
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
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
# Solve for free ligand from quasi-steady-state
L := 0.5 * ((Ltot - Rtot - KSS) + sqrt((Ltot - Rtot - KSS)^2 + 4 * KSS * Ltot))
# Complex and free receptor from mass balance
LR := Ltot - L
R := Rtot - LR
end
@dynamics begin
Ltot' = -Kel * L - Kint * LR
Rtot' = Ksyn - Kdeg * R - Kint * LR
end
@derived begin
Free_Ligand = @. L
Free_Receptor = @. R
Complex_Conc = @. LR
Total_Ligand = @. Ltot
Total_Receptor = @. Rtot
end
endPumasModel
Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKSS, tvKint
Random effects:
Covariates:
Dynamical system variables: Ltot, Rtot
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
Observed: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor
8.4 Relationship Between Full Model and QSS
The full model parameters can be recovered from QSS parameters if both \(K_D\) and \(k_{int}\) are known:
\[ K_D = K_{SS} - \frac{k_{int}}{k_{on}} \]
Or equivalently:
\[ k_{on} = \frac{k_{int}}{K_{SS} - K_D} \]
In practice, QSS is often used when \(k_{on}\) and \(k_{off}\) cannot be separately identified.
9 Michaelis-Menten (MM) Approximation
9.1 Mathematical Basis
The Michaelis-Menten approximation is the most simplified TMDD model. It applies when the ligand concentration is much greater than the total receptor concentration:
\[ L \gg R_{tot} \]
Under this assumption, the target is essentially fully saturated at therapeutic concentrations, and the target-mediated elimination reduces to a capacity-limited (saturable) process:
\[ \frac{dL}{dt} = -k_{el} \cdot L - \frac{V_{max} \cdot L}{K_M + L} \]
The MM model can be derived from different parent models, leading to different interpretations of \(K_M\):
1. From QSS (reversible binding): \[V_{max} = R_0 \cdot k_{int}, \quad K_M = K_{SS} = \frac{k_{off} + k_{int}}{k_{on}}\]
2. From Irreversible Binding (Gibiansky & Gibiansky, 2009): \[V_{max} = k_{syn}, \quad K_M = K_{IB} = \frac{k_{deg}}{k_{on}}\]
The irreversible binding derivation applies when \(k_{int} \gg k_{off}\) (complex degradation much faster than dissociation). In practice, both derivations yield similar MM kinetics, but the parameter relationships differ.
9.2 When is MM Appropriate?
The Michaelis-Menten approximation is valid when:
- Free receptor and complex concentrations are not measured
- Ligand concentrations are always above target concentrations
- Target-mediated clearance follows saturable kinetics
- Only ligand PK data are available
The MM approximation cannot describe receptor dynamics. If target engagement or receptor occupancy predictions are needed, use the QSS or full model.
9.3 MM Model in Pumas
tmdd_mm_model = @model begin
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvVmax ∈ RealDomain(lower = 0) # Maximum elimination rate
tvKM ∈ RealDomain(lower = 0) # Michaelis constant
end
@pre begin
Cl = tvCl
Vc = tvVc
Vmax = tvVmax
KM = tvKM
end
@dosecontrol begin
bioav = (Ligand = 1 / tvVc,)
end
@vars begin
Kel := Cl / Vc
end
@dynamics begin
Ligand' = -Kel * Ligand - Vmax * Ligand / (KM + Ligand)
end
@derived begin
Free_Ligand = @. Ligand
end
endPumasModel
Parameters: tvCl, tvVc, tvVmax, tvKM
Random effects:
Covariates:
Dynamical system variables: Ligand
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand
Observed: Free_Ligand
The MM model has only 4 parameters (\(CL\), \(V_c\), \(V_{max}\), \(K_M\)), making it the most identifiable but least mechanistic approximation. It is often a good starting point for model development.
9.4 Terminal Slope Limitations of MM
A critical insight from Peletier & Gabrielsson (2012) is that the Michaelis-Menten model can overestimate the terminal elimination rate compared to the full TMDD model.
9.4.1 Why This Matters
At low ligand concentrations (terminal phase), the elimination behavior differs:
Full TMDD Model (internalization-limited):
\[ \lambda_z^{TMDD} \approx k_{int} \]
When receptor turnover is slow relative to internalization, the terminal slope is governed by the internalization rate of the drug-receptor complex.
Michaelis-Menten Model:
\[ \lambda_z^{MM} = k_{el} + \frac{V_{max}}{K_M \cdot V_c} = k_{el} + \frac{k_{int} \cdot R_0}{K_M} \]
The MM model includes an additional term that can significantly increase the apparent terminal slope.
If you fit a Michaelis-Menten model to TMDD data and use the terminal slope to predict long-term drug exposure, you may underestimate the drug concentrations at late time points. This is because the MM approximation assumes \(L \gg R_{tot}\) throughout, which breaks down in the terminal phase.
9.4.2 When MM Terminal Slope is Accurate
The MM approximation for terminal slope is accurate when:
- \(k_{el} \gg k_{int}\) (linear elimination dominates)
- \(R_0 \ll K_M\) (low receptor relative to binding constant)
- Data only extend to the transition phase (phase C) and don’t capture the true terminal phase (phase D)
9.4.3 Quantifying the Difference
The discrepancy ratio is approximately:
\[ \frac{\lambda_z^{MM}}{\lambda_z^{TMDD}} \approx \frac{k_{el}}{k_{int}} + \frac{R_0}{K_M} \]
For a typical mAb with \(k_{el} \approx 0.002\) hr\(^{-1}\), \(k_{int} \approx 0.002\) hr\(^{-1}\), \(R_0 = 5\) nM, and \(K_M = 1\) nM:
\[ \text{Ratio} \approx 1 + 5 = 6 \]
This means the MM model might predict a terminal half-life 6 times shorter than reality!
When accurate prediction of terminal phase kinetics is important (e.g., for washout periods, immunogenicity assessment, or exposure at trough), prefer the QSS or full TMDD model over Michaelis-Menten.
10 Parameter Comparison Across Models
The following table shows how parameters relate across approximation levels:
| Full Model | QE Model | QSS Model | MM Model | Relationship |
|---|---|---|---|---|
| \(k_{on}\), \(k_{off}\) | \(K_D\) | \(K_{SS}\) | \(K_M\) | \(K_D = k_{off}/k_{on}\); \(K_{SS} = (k_{off}+k_{int})/k_{on}\) |
| \(k_{syn}\), \(k_{deg}\) | \(k_{syn}\), \(k_{deg}\) | \(k_{syn}\), \(k_{deg}\) | Implicit in \(V_{max}\) | \(R_0 = k_{syn}/k_{deg}\) |
| \(k_{int}\) | \(k_{int}\) | \(k_{int}\) | Implicit in \(V_{max}\) | \(V_{max} = R_0 \cdot k_{int}\) |
| \(CL\), \(V_c\) | \(CL\), \(V_c\) | \(CL\), \(V_c\) | \(CL\), \(V_c\) | Same |
11 Comparing Model Predictions
Let’s simulate all four models with equivalent parameters to visualize their predictions.
Show/Hide Code
# Define common parameters based on typical mAb values
R0 = 5.0 # Baseline receptor (nM)
Kdeg = 0.008 # Receptor degradation (hr⁻¹) ~0.2 day⁻¹
Ksyn = R0 * Kdeg # Receptor synthesis (nM/hr)
Kon = 0.33 # Association rate (nM⁻¹ hr⁻¹) ~8 day⁻¹
Koff = 0.33 # Dissociation rate (hr⁻¹) ~8 day⁻¹
Kint = 0.002 # Internalization rate (hr⁻¹) ~0.04 day⁻¹
KD = Koff / Kon # = 1.0 nM
KSS = (Koff + Kint) / Kon # ≈ 1.006 nM
Vmax = R0 * Kint # = 0.01 nM/hr
# Parameter sets for each model
param_full = (
tvCl = 0.006,
tvVc = 3.0,
tvKsyn = Ksyn,
tvKdeg = Kdeg,
tvKon = Kon,
tvKoff = Koff,
tvKint = Kint,
)
param_qe =
(tvCl = 0.006, tvVc = 3.0, tvKsyn = Ksyn, tvKdeg = Kdeg, tvKD = KD, tvKint = Kint)
param_qss =
(tvCl = 0.006, tvVc = 3.0, tvKsyn = Ksyn, tvKdeg = Kdeg, tvKSS = KSS, tvKint = Kint)
param_mm = (tvCl = 0.006, tvVc = 3.0, tvVmax = Vmax, tvKM = KSS)
# Create population with varying doses (wide range for clear separation)
doses = [5, 50, 500]
Random.seed!(42)
pop_full = [
Subject(
id = i,
events = DosageRegimen(dose, time = 0, cmt = 1),
observations = (
Free_Ligand = nothing,
Free_Receptor = nothing,
Complex_Conc = nothing,
),
) for (i, dose) in enumerate(doses)
]
pop_qe = [
Subject(
id = i,
events = DosageRegimen(dose, time = 0, cmt = 1),
observations = (
Free_Ligand = nothing,
Free_Receptor = nothing,
Complex_Conc = nothing,
),
) for (i, dose) in enumerate(doses)
]
pop_qss = [
Subject(
id = i,
events = DosageRegimen(dose, time = 0, cmt = 1),
observations = (
Free_Ligand = nothing,
Free_Receptor = nothing,
Complex_Conc = nothing,
),
) for (i, dose) in enumerate(doses)
]
pop_mm = [
Subject(
id = i,
events = DosageRegimen(dose, time = 0, cmt = 1),
observations = (Free_Ligand = nothing,),
) for (i, dose) in enumerate(doses)
]
# Simulate all models (extended time for mAb kinetics)
obstimes = 0.1:2:672
sim_full = simobs(tmdd_full_model, pop_full, param_full, obstimes = obstimes)
sim_qe = simobs(tmdd_qe_model, pop_qe, param_qe, obstimes = obstimes)
sim_qss = simobs(tmdd_qss_model, pop_qss, param_qss, obstimes = obstimes)
sim_mm = simobs(tmdd_mm_model, pop_mm, param_mm, obstimes = obstimes)Simulated population (Vector{<:Subject})
Simulated subjects: 3
Simulated variables: Free_Ligand
11.1 Free Ligand Comparison
Show/Hide Code
# Convert all simulations to DataFrames and combine
df_full = DataFrame(sim_full)
df_full.Model .= "Full"
df_qe = DataFrame(sim_qe)
df_qe.Model .= "QE"
df_qss = DataFrame(sim_qss)
df_qss.Model .= "QSS"
df_mm = DataFrame(sim_mm)
df_mm.Model .= "MM"
# Create dose labels
dose_label_df = DataFrame(id = string.(1:3), Dose = string.(doses) .* " units")
df_full = innerjoin(df_full, dose_label_df, on = :id)
df_qe = innerjoin(df_qe, dose_label_df, on = :id)
df_qss = innerjoin(df_qss, dose_label_df, on = :id)
df_mm = innerjoin(df_mm, dose_label_df, on = :id)
# Combine all models
df_models = vcat(
select(df_full, :time, :Free_Ligand, :Model, :Dose),
select(df_qe, :time, :Free_Ligand, :Model, :Dose),
select(df_qss, :time, :Free_Ligand, :Model, :Dose),
select(df_mm, :time, :Free_Ligand, :Model, :Dose),
)
plt =
data(df_models) *
mapping(
:time => "Time (hr)",
:Free_Ligand => "Free Ligand Concentration (nM)";
color = :Model,
linestyle = :Model,
layout = :Dose,
) *
visual(Lines; linewidth = 2)
draw(plt; axis = (yscale = Makie.pseudolog10,))11.2 Receptor and Complex Dynamics
The MM model cannot predict receptor or complex dynamics. Let’s compare the other three models:
Show/Hide Code
# Use middle dose for comparison
dose_idx = 2
# Prepare data for models that have receptor/complex
df_receptor = vcat(
select(
filter(row -> row.id == string(dose_idx), df_full),
:time,
:Free_Receptor,
:Complex_Conc,
:Model,
),
select(
filter(row -> row.id == string(dose_idx), df_qe),
:time,
:Free_Receptor,
:Complex_Conc,
:Model,
),
select(
filter(row -> row.id == string(dose_idx), df_qss),
:time,
:Free_Receptor,
:Complex_Conc,
:Model,
),
)
# Stack for faceted plot
df_long = stack(
df_receptor,
[:Free_Receptor, :Complex_Conc],
variable_name = :Species,
value_name = :Concentration,
)
df_long.Species =
replace.(
string.(df_long.Species),
"Free_Receptor" => "Free Receptor",
"Complex_Conc" => "Complex",
)
plt =
data(df_long) *
mapping(
:time => "Time (hr)",
:Concentration => "Concentration (nM)";
color = :Model,
linestyle = :Model,
layout = :Species,
) *
visual(Lines; linewidth = 2)
draw(plt; figure = (title = "Dose = $(doses[dose_idx]) units",))12 Model Selection Decision Framework
Choosing the right TMDD approximation depends on several factors. Here is a practical decision framework:
12.1 Step 1: What Data Are Available?
| Data Available | Recommended Starting Model |
|---|---|
| Free ligand only | Michaelis-Menten |
| Free ligand + total receptor | QSS |
| Free ligand + total ligand | QSS |
| Free ligand + free receptor + complex | Full TMDD |
12.2 Step 2: Consider Scientific Context
- Do you need receptor dynamics? Use QSS or Full model
- Is internalization significant? Use QSS over QE
- Is target concentration known to be low relative to drug? MM may suffice
- Do you need mechanistic interpretation? Use Full or QSS model
12.3 Step 3: Assess Model Fit and Identifiability
Start with the simplest appropriate model (often MM) and increase complexity if:
- The model fails to capture key PK features
- Additional data become available
- Mechanistic predictions are needed
Start simple, add complexity as needed.
Fit Michaelis-Menten model first
If MM fails, try QSS
If QSS fails and you have receptor data, try Full TMDD
Compare models using information criteria (AIC, BIC)
:::
13 Identifiability Considerations
Following Gibiansky et al. (2008), here are key identifiability insights:
13.1 With Free Ligand Data Only
- MM model: All 4 parameters identifiable
- QSS model: \(CL\), \(V_c\) identifiable; \(k_{syn}\), \(k_{deg}\), \(K_{SS}\), \(k_{int}\) may require constraints
- Full model: Generally not identifiable
13.2 With Total Ligand and Total Receptor Data
- QSS model: All parameters typically identifiable
- Full model: Still may have issues with \(k_{on}\) and \(k_{off}\) separately
13.3 Practical Strategies
- Fix known parameters: Use literature values for \(k_{deg}\), \(K_D\) from in vitro studies
- Reduce model: If parameters are not identifiable, use a simpler approximation
- Collect more data: Measure receptor concentrations when possible
- Profile likelihood: Assess practical identifiability of estimated parameters
14 Adding Peripheral Compartment (Two-Compartment Models)
All TMDD models can be extended to include peripheral distribution. Here is the two-compartment QSS model:
tmdd_qss_2cmt_model = @model begin
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvQ ∈ RealDomain(lower = 0) # Intercompartmental clearance
tvVp ∈ RealDomain(lower = 0) # Peripheral volume
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKSS ∈ RealDomain(lower = 0)
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
Q = tvQ
Vp = tvVp
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
K12 := Q / Vc
K21 := Q / Vp
R0 := Ksyn / Kdeg
# Solve for free ligand from quasi-steady-state
L := 0.5 * ((Ltot - Rtot - KSS) + sqrt((Ltot - Rtot - KSS)^2 + 4 * KSS * Ltot))
# Complex and free receptor from mass balance
LR := Ltot - L
R := Rtot - LR
end
@dynamics begin
Ltot' = -Kel * L - Kint * LR - K12 * L + K21 * (Peripheral / Vc)
Peripheral' = K12 * L * Vc - K21 * Peripheral
Rtot' = Ksyn - Kdeg * R - Kint * LR
end
@derived begin
Free_Ligand = @. L
Free_Receptor = @. R
Complex_Conc = @. LR
end
endPumasModel
Parameters: tvCl, tvVc, tvQ, tvVp, tvKsyn, tvKdeg, tvKSS, tvKint
Random effects:
Covariates:
Dynamical system variables: Ltot, Peripheral, Rtot
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc
Observed: Free_Ligand, Free_Receptor, Complex_Conc
15 Summary
In this tutorial, we covered:
- Why simplify: Identifiability, data limitations, numerical stability
- QE approximation: Assumes binding at equilibrium; uses \(K_D\)
- QSS approximation: More general; uses \(K_{SS} = (k_{off} + k_{int})/k_{on}\)
- MM approximation: Simplest; assumes \(L \gg R_{tot}\); uses \(V_{max}\) and \(K_M\)
- Model selection: Start simple, add complexity based on data and needs
- Identifiability: Key consideration for parameter estimation
In Part 4, we will explore special cases and extensions including constant receptor models, irreversible binding, and multi-target scenarios.