using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using Random
TMDD Part 4: Special Cases and Extensions
1 Introduction
This is Part 4 of our comprehensive TMDD tutorial series. Here we explore specialized TMDD scenarios that extend beyond the standard model, including constant receptor models, irreversible binding, drug-drug interactions at the target level, and FcRn-mediated recycling.
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 (this tutorial) - 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:
- Implement the constant receptor (\(R_{tot}\)) TMDD model
- Model irreversible (covalent) binding scenarios
- Handle two drugs competing for the same target
- Incorporate FcRn-mediated recycling for mAbs
- Model subcutaneous administration with absorption-site binding
- Understand when each special case applies
3 Libraries
4 Constant Total Receptor Model
4.1 When to Use
The constant \(R_{tot}\) model (Wagner assumption) applies when:
- Target turnover is rapid compared to drug disposition
- Receptor is continuously replenished from a large pool
- Receptor synthesis rate is high relative to drug binding
- Critically: \(k_{int} \approx k_{deg}\) (the complex and free receptor degrade at similar rates)
As noted in the literature (Gibiansky et al., 2008; Wagner, 1976), the constant \(R_{tot}\) assumption requires that \(k_{int} \approx k_{deg}\). When this condition holds, the total receptor (\(R + LR\)) remains constant because the degradation rate of the complex equals that of the free receptor.
Under these conditions, the total receptor concentration remains approximately constant:
\[ R_{tot} = R + LR \approx \text{constant} \]
4.2 Mathematical Formulation
Starting from the standard TMDD model and assuming \(R_{tot}\) is constant:
\[ \frac{dL}{dt} = -k_{el} \cdot L - k_{on} \cdot L \cdot (R_{tot} - LR) + k_{off} \cdot LR \]
\[ \frac{dLR}{dt} = k_{on} \cdot L \cdot (R_{tot} - LR) - k_{off} \cdot LR - k_{int} \cdot LR \]
Or equivalently:
\[ \frac{dLR}{dt} = k_{on} \cdot L \cdot R_{tot} - (k_{on} \cdot L + k_{off} + k_{int}) \cdot LR \]
4.3 Pumas Implementation
tmdd_constant_rtot = @model begin
@metadata begin
desc = "TMDD with Constant Total Receptor"
timeu = u"hr"
end
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvRtot ∈ RealDomain(lower = 0) # Constant total receptor
tvKon ∈ RealDomain(lower = 0)
tvKoff ∈ RealDomain(lower = 0)
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
Rtot = tvRtot # Constant
Kon = tvKon
Koff = tvKoff
Kint = tvKint
end
@dosecontrol begin
bioav = (Ligand = 1 / tvVc,)
end
@vars begin
Kel := Cl / Vc
KD := Koff / Kon
# Free receptor from mass balance (Rtot is constant)
R := Rtot - Complex
end
@dynamics begin
# Free ligand
Ligand' = -Kel * Ligand - Kon * Ligand * R + Koff * Complex
# Complex (no receptor dynamics needed since Rtot is constant)
Complex' = Kon * Ligand * R - Koff * Complex - Kint * Complex
end
@derived begin
Free_Ligand = @. Ligand
Free_Receptor = @. R
Complex_Conc = @. Complex
Target_Occupancy = @. Complex / Rtot
end
end┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Free_Receptor` is defined in the `@derived` block using `=` and hence `Free_Receptor` is not used for model fitting but only returned when simulating: │ If `Free_Receptor` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Receptor` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Receptor` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Complex_Conc` is defined in the `@derived` block using `=` and hence `Complex_Conc` is not used for model fitting but only returned when simulating: │ If `Complex_Conc` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Complex_Conc` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Complex_Conc` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Target_Occupancy` is defined in the `@derived` block using `=` and hence `Target_Occupancy` is not used for model fitting but only returned when simulating: │ If `Target_Occupancy` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Target_Occupancy` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Target_Occupancy` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
PumasModel
Parameters: tvCl, tvVc, tvRtot, tvKon, tvKoff, tvKint
Random effects:
Covariates:
Dynamical system variables: Ligand, Complex
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc, Target_Occupancy
Observed: Free_Ligand, Free_Receptor, Complex_Conc, Target_Occupancy
4.4 Simulation Example
Show/Hide Code
# Parameters for constant Rtot model
# Note: kint ≈ kdeg for the constant Rtot assumption to hold
param_const_rtot = (
tvCl = 0.006,
tvVc = 3.0,
tvRtot = 5.0, # Constant at 5 nM (similar to R0 in other tutorials)
tvKon = 0.33,
tvKoff = 0.33,
tvKint = 0.008, # Set equal to kdeg (0.008 hr⁻¹) for constant Rtot
)
doses = [5.0, 50.0, 500.0] # nmol - wide range for clear visualization
Random.seed!(42)
pop_const = [
Subject(
id = i,
events = DosageRegimen(dose, time = 0, cmt = 1),
observations = (
Free_Ligand = nothing,
Free_Receptor = nothing,
Target_Occupancy = nothing,
),
) for (i, dose) in enumerate(doses)
]
obstimes = 0.1:2:672 # Extended time for mAb kinetics
sim_const = simobs(tmdd_constant_rtot, pop_const, param_const_rtot, obstimes = obstimes)Simulated population (Vector{<:Subject})
Simulated subjects: 3
Simulated variables: Free_Ligand, Free_Receptor, Complex_Conc, Target_Occupancy
Show/Hide Code
# Convert simulation to DataFrame
df_const = DataFrame(sim_const)
dose_label_df = DataFrame(id = string.(1:3), Dose = ["5 nmol", "50 nmol", "500 nmol"])
df_const = innerjoin(df_const, dose_label_df, on = :id)
# Stack for faceted plot
df_long = stack(
df_const,
[:Free_Receptor, :Target_Occupancy],
variable_name = :Metric,
value_name = :Value,
)
df_long.Metric =
replace.(
string.(df_long.Metric),
"Free_Receptor" => "Free Receptor (nM)",
"Target_Occupancy" => "Target Occupancy",
)
plt =
data(df_long) *
mapping(:time => "Time (hr)", :Value; color = :Dose, layout = :Metric) *
visual(Lines; linewidth = 2)
draw(plt; figure = (title = "Constant Rtot Model",))With the constant \(R_{tot}\) model, the free receptor returns to exactly \(R_{tot}\) after drug washout, regardless of dose. This contrasts with the standard model where receptor can be temporarily depleted.
5 Irreversible Binding Model
5.1 When to Use
Irreversible binding (\(k_{off} = 0\)) applies for:
- Covalent inhibitors that form permanent bonds
- Very tight binders where \(k_{off}\) is negligibly small (\(k_{off} \ll k_{int}, k_{deg}\))
- Suicide substrates that inactivate the target
5.2 Mathematical Formulation
With \(k_{off} = 0\):
\[ \frac{dL}{dt} = -k_{el} \cdot L - k_{on} \cdot L \cdot R \]
\[ \frac{dR}{dt} = k_{syn} - k_{deg} \cdot R - k_{on} \cdot L \cdot R \]
\[ \frac{dLR}{dt} = k_{on} \cdot L \cdot R - k_{int} \cdot LR \]
Notice: the complex cannot dissociate, only be internalized/degraded.
For irreversible binding models, Gibiansky & Gibiansky (2009) introduced the irreversible binding constant:
\[ K_{IB} = \frac{k_{deg}}{k_{on}} \]
This parameter plays an analogous role to \(K_D\) in reversible binding. When the irreversible binding model is further reduced to Michaelis-Menten, \(K_M = K_{IB}\) and \(V_{max} = k_{syn}\).
The relationship \(K_{IB} = k_{deg}/k_{on}\) differs from \(K_D = k_{off}/k_{on}\) because in irreversible binding, the receptor turnover (\(k_{deg}\)) determines the effective binding capacity rather than dissociation.
5.3 Pumas Implementation
tmdd_irreversible = @model begin
@metadata begin
desc = "TMDD with Irreversible Binding (Koff = 0)"
timeu = u"hr"
end
@param begin
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKon ∈ RealDomain(lower = 0) # Association only
tvKint ∈ RealDomain(lower = 0)
end
@pre begin
Cl = tvCl
Vc = tvVc
Ksyn = tvKsyn
Kdeg = tvKdeg
Kon = tvKon
Kint = tvKint
end
@dosecontrol begin
bioav = (Ligand = 1 / tvVc,)
end
@init begin
Receptor = Ksyn / Kdeg
end
@vars begin
Kel := Cl / Vc
R0 := Ksyn / Kdeg
end
@dynamics begin
# No Koff term - binding is irreversible
Ligand' = -Kel * Ligand - Kon * Ligand * Receptor
Receptor' = Ksyn - Kdeg * Receptor - Kon * Ligand * Receptor
# Complex only cleared by internalization
Complex' = Kon * Ligand * Receptor - Kint * Complex
end
@derived begin
Free_Ligand = @. Ligand
Free_Receptor = @. Receptor
Complex_Conc = @. Complex
end
end┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Free_Receptor` is defined in the `@derived` block using `=` and hence `Free_Receptor` is not used for model fitting but only returned when simulating: │ If `Free_Receptor` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Receptor` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Receptor` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Complex_Conc` is defined in the `@derived` block using `=` and hence `Complex_Conc` is not used for model fitting but only returned when simulating: │ If `Complex_Conc` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Complex_Conc` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Complex_Conc` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
PumasModel
Parameters: tvCl, tvVc, tvKsyn, tvKdeg, tvKon, tvKint
Random effects:
Covariates:
Dynamical system variables: Ligand, Receptor, Complex
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Free_Receptor, Complex_Conc
Observed: Free_Ligand, Free_Receptor, Complex_Conc
5.4 Comparison: Reversible vs Irreversible
Show/Hide Code
# Define comparison models using consistent Gibiansky-based parameters
param_reversible = (
tvCl = 0.006,
tvVc = 3.0,
tvKsyn = 0.04,
tvKdeg = 0.008,
tvKon = 0.33,
tvKoff = 0.33, # Reversible (KD = 1 nM)
tvKint = 0.002,
)
param_irreversible =
(tvCl = 0.006, tvVc = 3.0, tvKsyn = 0.04, tvKdeg = 0.008, tvKon = 0.33, tvKint = 0.002)
# Standard TMDD model for comparison
tmdd_reversible = @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
@dynamics begin
Ligand' = -(Cl / Vc) * 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
end
end
single_pop = [
Subject(
id = 1,
events = DosageRegimen(50, time = 0, cmt = 1), # 50 nmol dose
observations = (Free_Ligand = nothing, Free_Receptor = nothing),
),
]
sim_rev = simobs(tmdd_reversible, single_pop, param_reversible, obstimes = obstimes)
sim_irrev = simobs(tmdd_irreversible, single_pop, param_irreversible, obstimes = obstimes)
# Convert to DataFrames
df_rev = DataFrame(sim_rev)
df_rev.Binding .= "Reversible"
df_irrev = DataFrame(sim_irrev)
df_irrev.Binding .= "Irreversible"
df_binding = vcat(
select(df_rev, :time, :Free_Ligand, :Free_Receptor, :Binding),
select(df_irrev, :time, :Free_Ligand, :Free_Receptor, :Binding),
)
# Stack for faceted plot
df_long = stack(
df_binding,
[:Free_Ligand, :Free_Receptor],
variable_name = :Species,
value_name = :Concentration,
)
df_long.Species =
replace.(
string.(df_long.Species),
"Free_Ligand" => "Free Ligand",
"Free_Receptor" => "Free Receptor",
)
plt =
data(df_long) *
mapping(
:time => "Time (hr)",
:Concentration => "Concentration (nM)";
color = :Binding,
linestyle = :Binding,
layout = :Species,
) *
visual(Lines; linewidth = 2)
draw(plt)┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Free_Receptor` is defined in the `@derived` block using `=` and hence `Free_Receptor` is not used for model fitting but only returned when simulating: │ If `Free_Receptor` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Receptor` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Receptor` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
With irreversible binding:
- Free ligand clears faster (no release from complex)
- Receptor suppression is more pronounced
- Recovery depends entirely on new receptor synthesis
6 Two Drugs Competing for the Same Target
6.1 When to Use
Competition models apply for:
- Drug-drug interactions at the target level
- Biosimilar comparisons with the originator
- Combination therapies targeting the same receptor
- Displacement studies
6.2 Mathematical Formulation
The competitive binding equations, originally described by Gaddum (1926), extend the TMDD framework for two ligands:
For two ligands (\(L_1\) and \(L_2\)) competing for the same receptor:
\[ \frac{dL_1}{dt} = -k_{el,1} \cdot L_1 - k_{on,1} \cdot L_1 \cdot R + k_{off,1} \cdot L_1R \]
\[ \frac{dL_2}{dt} = -k_{el,2} \cdot L_2 - k_{on,2} \cdot L_2 \cdot R + k_{off,2} \cdot L_2R \]
\[ \frac{dR}{dt} = k_{syn} - k_{deg} \cdot R - k_{on,1} \cdot L_1 \cdot R + k_{off,1} \cdot L_1R - k_{on,2} \cdot L_2 \cdot R + k_{off,2} \cdot L_2R \]
\[ \frac{dL_1R}{dt} = k_{on,1} \cdot L_1 \cdot R - k_{off,1} \cdot L_1R - k_{int,1} \cdot L_1R \]
\[ \frac{dL_2R}{dt} = k_{on,2} \cdot L_2 \cdot R - k_{off,2} \cdot L_2R - k_{int,2} \cdot L_2R \]
6.3 Pumas Implementation
tmdd_competition = @model begin
@metadata begin
desc = "TMDD with Two Competing Ligands"
timeu = u"hr"
end
@param begin
# Drug 1 parameters
tvCl1 ∈ RealDomain(lower = 0)
tvVc1 ∈ RealDomain(lower = 0)
tvKon1 ∈ RealDomain(lower = 0)
tvKoff1 ∈ RealDomain(lower = 0)
tvKint1 ∈ RealDomain(lower = 0)
# Drug 2 parameters
tvCl2 ∈ RealDomain(lower = 0)
tvVc2 ∈ RealDomain(lower = 0)
tvKon2 ∈ RealDomain(lower = 0)
tvKoff2 ∈ RealDomain(lower = 0)
tvKint2 ∈ RealDomain(lower = 0)
# Receptor parameters
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
end
@pre begin
Cl1 = tvCl1
Vc1 = tvVc1
Kon1 = tvKon1
Koff1 = tvKoff1
Kint1 = tvKint1
Cl2 = tvCl2
Vc2 = tvVc2
Kon2 = tvKon2
Koff2 = tvKoff2
Kint2 = tvKint2
Ksyn = tvKsyn
Kdeg = tvKdeg
end
@dosecontrol begin
bioav = (Ligand1 = 1 / tvVc1, Ligand2 = 1 / tvVc2)
end
@init begin
Receptor = Ksyn / Kdeg
end
@vars begin
Kel1 := Cl1 / Vc1
Kel2 := Cl2 / Vc2
R0 := Ksyn / Kdeg
Rtot := Receptor + Complex1 + Complex2
end
@dynamics begin
# Drug 1
Ligand1' = -Kel1 * Ligand1 - Kon1 * Ligand1 * Receptor + Koff1 * Complex1
# Drug 2
Ligand2' = -Kel2 * Ligand2 - Kon2 * Ligand2 * Receptor + Koff2 * Complex2
# Shared receptor - both drugs compete
Receptor' =
Ksyn - Kdeg * Receptor - Kon1 * Ligand1 * Receptor + Koff1 * Complex1 -
Kon2 * Ligand2 * Receptor + Koff2 * Complex2
# Complex 1 (Drug 1-Receptor)
Complex1' = Kon1 * Ligand1 * Receptor - Koff1 * Complex1 - Kint1 * Complex1
# Complex 2 (Drug 2-Receptor)
Complex2' = Kon2 * Ligand2 * Receptor - Koff2 * Complex2 - Kint2 * Complex2
end
@derived begin
Free_Ligand1 = @. Ligand1
Free_Ligand2 = @. Ligand2
Free_Receptor = @. Receptor
Occupancy_Drug1 = @. Complex1 / Rtot
Occupancy_Drug2 = @. Complex2 / Rtot
Total_Occupancy = @. (Complex1 + Complex2) / Rtot
end
end┌ Warning: Variable `Free_Ligand1` is defined in the `@derived` block using `=` and hence `Free_Ligand1` is not used for model fitting but only returned when simulating: │ If `Free_Ligand1` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand1` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand1` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Free_Ligand2` is defined in the `@derived` block using `=` and hence `Free_Ligand2` is not used for model fitting but only returned when simulating: │ If `Free_Ligand2` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand2` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand2` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Free_Receptor` is defined in the `@derived` block using `=` and hence `Free_Receptor` is not used for model fitting but only returned when simulating: │ If `Free_Receptor` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Receptor` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Receptor` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Occupancy_Drug1` is defined in the `@derived` block using `=` and hence `Occupancy_Drug1` is not used for model fitting but only returned when simulating: │ If `Occupancy_Drug1` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Occupancy_Drug1` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Occupancy_Drug1` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Occupancy_Drug2` is defined in the `@derived` block using `=` and hence `Occupancy_Drug2` is not used for model fitting but only returned when simulating: │ If `Occupancy_Drug2` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Occupancy_Drug2` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Occupancy_Drug2` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Total_Occupancy` is defined in the `@derived` block using `=` and hence `Total_Occupancy` is not used for model fitting but only returned when simulating: │ If `Total_Occupancy` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Total_Occupancy` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Total_Occupancy` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
PumasModel
Parameters: tvCl1, tvVc1, tvKon1, tvKoff1, tvKint1, tvCl2, tvVc2, tvKon2, tvKoff2, tvKint2, tvKsyn, tvKdeg
Random effects:
Covariates:
Dynamical system variables: Ligand1, Ligand2, Receptor, Complex1, Complex2
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand1, Free_Ligand2, Free_Receptor, Occupancy_Drug1, Occupancy_Drug2, Total_Occupancy
Observed: Free_Ligand1, Free_Ligand2, Free_Receptor, Occupancy_Drug1, Occupancy_Drug2, Total_Occupancy
6.4 Competition Simulation Example
Show/Hide Code
# Parameters: Drug 1 has higher affinity (lower KD)
# Using literature-based parameters
param_competition = (
tvCl1 = 0.006,
tvVc1 = 3.0,
tvKon1 = 0.5,
tvKoff1 = 0.1,
tvKint1 = 0.002, # High affinity (KD=0.2 nM)
tvCl2 = 0.006,
tvVc2 = 3.0,
tvKon2 = 0.33,
tvKoff2 = 0.66,
tvKint2 = 0.002, # Lower affinity (KD=2 nM)
tvKsyn = 0.04,
tvKdeg = 0.008,
)
# Scenario: Drug 2 given first, Drug 1 added later
regimen_comp = DosageRegimen([
DosageRegimen(50, time = 0, cmt = 2), # Drug 2 at t=0 (50 nmol)
DosageRegimen(50, time = 168, cmt = 1), # Drug 1 at t=168 hr (1 week later)
])
pop_comp = [
Subject(
id = 1,
events = regimen_comp,
observations = (
Free_Ligand1 = nothing,
Free_Ligand2 = nothing,
Free_Receptor = nothing,
Occupancy_Drug1 = nothing,
Occupancy_Drug2 = nothing,
Total_Occupancy = nothing,
),
),
]
obstimes_comp = 0.1:2:672 # 4 weeks
sim_comp = simobs(tmdd_competition, pop_comp, param_competition, obstimes = obstimes_comp)
# Convert to DataFrame
df_comp = DataFrame(sim_comp)
# Panel 1: Free Drug Concentrations
df_drugs = stack(
df_comp,
[:Free_Ligand1, :Free_Ligand2],
variable_name = :Drug,
value_name = :Concentration,
)
df_drugs.Drug =
replace.(
string.(df_drugs.Drug),
"Free_Ligand1" => "Drug 1 (high affinity)",
"Free_Ligand2" => "Drug 2 (lower affinity)",
)
plt1 =
data(df_drugs) *
mapping(:time => "Time (hr)", :Concentration => "Concentration (nM)"; color = :Drug) *
visual(Lines; linewidth = 2)
# Panel 2: Free Receptor
plt2 =
data(df_comp) *
mapping(:time => "Time (hr)", :Free_Receptor => "Free Receptor (nM)") *
visual(Lines; linewidth = 2, color = :green)
# Panel 3: Target Occupancy by Drug
df_occ = stack(
df_comp,
[:Occupancy_Drug1, :Occupancy_Drug2],
variable_name = :Drug,
value_name = :Occupancy,
)
df_occ.Drug =
replace.(
string.(df_occ.Drug),
"Occupancy_Drug1" => "Drug 1",
"Occupancy_Drug2" => "Drug 2",
)
plt3 =
data(df_occ) *
mapping(
:time => "Time (hr)",
:Occupancy => "Target Occupancy (fraction)";
color = :Drug,
) *
visual(Lines; linewidth = 2)
# Panel 4: Total Occupancy
plt4 =
data(df_comp) *
mapping(:time => "Time (hr)", :Total_Occupancy => "Total Occupancy") *
visual(Lines; linewidth = 2, color = :purple)
fig = Figure(size = (900, 600))
draw!(
fig[1, 1],
plt1;
axis = (yscale = Makie.pseudolog10, title = "Free Drug Concentrations"),
)
draw!(fig[1, 2], plt2; axis = (title = "Free Receptor",))
draw!(fig[2, 1], plt3; axis = (title = "Target Occupancy by Each Drug",))
draw!(fig[2, 2], plt4; axis = (title = "Total Target Occupancy",))
# Add vertical lines for Drug 1 addition (at 168 hr = 1 week)
for i = 1:2, j = 1:2
vlines!(fig[i, j][1, 1], [168], color = :gray, linestyle = :dash)
end
figWhen Drug 1 (higher affinity) is added at 168 hours (1 week):
- Drug 2 occupancy decreases as it is displaced
- Drug 1 rapidly captures available receptor
- Total occupancy remains high due to the higher-affinity competitor
7 FcRn-Mediated Recycling
7.1 Why FcRn Matters
Monoclonal antibodies have unusually long half-lives (weeks) due to the neonatal Fc receptor (FcRn) salvage pathway:
- Antibodies are taken up by endothelial cells via pinocytosis
- In the acidic endosome, IgG binds to FcRn
- FcRn-bound IgG is recycled to the cell surface and released
- Unbound IgG is degraded in lysosomes
This recycling reduces the effective clearance of IgG antibodies.
7.2 Model Structure
The FcRn model adds a recycling compartment:
tmdd_fcrn = @model begin
@metadata begin
desc = "TMDD with FcRn-Mediated Recycling"
timeu = u"hr"
end
@param begin
# PK parameters
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
tvQ ∈ RealDomain(lower = 0)
tvVp ∈ RealDomain(lower = 0)
# TMDD parameters
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKSS ∈ RealDomain(lower = 0)
tvKint ∈ RealDomain(lower = 0)
# FcRn recycling parameters
tvKup ∈ RealDomain(lower = 0) # Pinocytosis uptake rate
tvFR ∈ RealDomain(lower = 0, upper = 1) # Fraction recycled
tvKrec ∈ RealDomain(lower = 0) # Recycling rate
end
@pre begin
Cl = tvCl
Vc = tvVc
Q = tvQ
Vp = tvVp
Ksyn = tvKsyn
Kdeg = tvKdeg
KSS = tvKSS
Kint = tvKint
Kup = tvKup
FR = tvFR
Krec = tvKrec
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
# QSS binding
L := 0.5 * ((Ltot - Rtot - KSS) + sqrt((Ltot - Rtot - KSS)^2 + 4 * KSS * Ltot))
LR := Ltot - L
R := Rtot - LR
# Degradation rate for non-recycled fraction
Kdeg_endo := Krec * (1 - FR) / FR
end
@dynamics begin
# Central total ligand with FcRn recycling
Ltot' =
-Kel * L - Kint * LR - K12 * L + K21 * (Peripheral / Vc) - Kup * L +
Krec * Endosomal
# Peripheral
Peripheral' = K12 * L * Vc - K21 * Peripheral
# Receptor dynamics
Rtot' = Ksyn - Kdeg * R - Kint * LR
# Endosomal compartment (FcRn-bound)
Endosomal' = Kup * L * FR - Krec * Endosomal
end
@derived begin
Free_Ligand = @. L
Endosomal_Conc = @. Endosomal
Total_Receptor = @. Rtot
end
end┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Endosomal_Conc` is defined in the `@derived` block using `=` and hence `Endosomal_Conc` is not used for model fitting but only returned when simulating: │ If `Endosomal_Conc` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Endosomal_Conc` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Endosomal_Conc` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Total_Receptor` is defined in the `@derived` block using `=` and hence `Total_Receptor` is not used for model fitting but only returned when simulating: │ If `Total_Receptor` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Total_Receptor` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Total_Receptor` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
PumasModel
Parameters: tvCl, tvVc, tvQ, tvVp, tvKsyn, tvKdeg, tvKSS, tvKint, tvKup, tvFR, tvKrec
Random effects:
Covariates:
Dynamical system variables: Ltot, Peripheral, Rtot, Endosomal
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Endosomal_Conc, Total_Receptor
Observed: Free_Ligand, Endosomal_Conc, Total_Receptor
7.3 Impact of FcRn Recycling
Show/Hide Code
# Model without FcRn recycling (standard TMDD)
tmdd_no_fcrn = @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)
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
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
Free_Ligand = @. L
end
end
# Parameters using literature-based values
param_no_fcrn =
(tvCl = 0.006, tvVc = 3.0, tvKsyn = 0.04, tvKdeg = 0.008, tvKSS = 1.0, tvKint = 0.002)
param_fcrn = (
tvCl = 0.006,
tvVc = 3.0,
tvQ = 0.02,
tvVp = 3.0,
tvKsyn = 0.04,
tvKdeg = 0.008,
tvKSS = 1.0,
tvKint = 0.002,
tvKup = 0.01, # Uptake rate
tvFR = 0.7, # 70% recycled
tvKrec = 0.02, # Recycling rate
)
pop_fcrn = [
Subject(
id = 1,
events = DosageRegimen(50, time = 0, cmt = 1),
observations = (Free_Ligand = nothing,),
),
]
obstimes_fcrn = 0.1:2:1008 # 6 weeks
sim_no_fcrn = simobs(tmdd_no_fcrn, pop_fcrn, param_no_fcrn, obstimes = obstimes_fcrn)
sim_fcrn = simobs(tmdd_fcrn, pop_fcrn, param_fcrn, obstimes = obstimes_fcrn)
# Convert to DataFrames
df_no_fcrn = DataFrame(sim_no_fcrn)
df_no_fcrn.Model .= "Without FcRn recycling"
df_fcrn = DataFrame(sim_fcrn)
df_fcrn.Model .= "With FcRn recycling (70%)"
df_fcrn_combined = vcat(
select(df_no_fcrn, :time, :Free_Ligand, :Model),
select(df_fcrn, :time, :Free_Ligand, :Model),
)
plt =
data(df_fcrn_combined) *
mapping(:time => "Time (hr)", :Free_Ligand => "Free Ligand (nM)"; color = :Model) *
visual(Lines; linewidth = 2)
draw(plt; axis = (yscale = Makie.pseudolog10, title = "Effect of FcRn Recycling on mAb PK"))┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
FcRn recycling:
- Extends terminal half-life significantly
- Reduces effective clearance
- Is saturable at very high doses (not modeled here)
- Can be engineered via Fc mutations (half-life extension technologies)
8 Subcutaneous Administration with Depot Binding
8.1 When to Use
For subcutaneous (SC) administration, target binding may occur at the injection site, affecting both bioavailability and absorption kinetics:
tmdd_sc = @model begin
@metadata begin
desc = "TMDD with SC Administration and Depot Binding"
timeu = u"hr"
end
@param begin
# Absorption parameters
tvKa ∈ RealDomain(lower = 0) # Absorption rate
tvF ∈ RealDomain(lower = 0, upper = 1) # Bioavailability
# PK parameters
tvCl ∈ RealDomain(lower = 0)
tvVc ∈ RealDomain(lower = 0)
# Central TMDD
tvKsyn ∈ RealDomain(lower = 0)
tvKdeg ∈ RealDomain(lower = 0)
tvKSS ∈ RealDomain(lower = 0)
tvKint ∈ RealDomain(lower = 0)
# Depot binding
tvRdepot ∈ RealDomain(lower = 0) # Target concentration at SC site
tvKSS_depot ∈ RealDomain(lower = 0) # Binding affinity at depot
end
@pre begin
Ka = tvKa
F = tvF
Cl = tvCl
Vc = tvVc
Ksyn = tvKsyn
Kdeg = tvKdeg
KSS = tvKSS
Kint = tvKint
Rdepot = tvRdepot
KSS_depot = tvKSS_depot
end
@dosecontrol begin
bioav = (Depot = tvF,)
end
@init begin
Rtot_central = Ksyn / Kdeg
end
@vars begin
Kel := Cl / Vc
R0 := Ksyn / Kdeg
# Central compartment QSS binding
L :=
0.5 * (
(Ltot - Rtot_central - KSS) +
sqrt((Ltot - Rtot_central - KSS)^2 + 4 * KSS * Ltot)
)
LR := Ltot - L
R := Rtot_central - LR
# Depot binding reduces free drug available for absorption
Depot_free := Depot / (1 + Rdepot / KSS_depot)
end
@dynamics begin
# SC depot with binding
Depot' = -Ka * Depot_free
# Central total ligand
Ltot' = Ka * Depot_free / Vc - Kel * L - Kint * LR
# Receptor turnover
Rtot_central' = Ksyn - Kdeg * R - Kint * LR
end
@derived begin
Free_Ligand = @. L
Depot_Amount = @. Depot
end
end┌ Warning: Variable `Free_Ligand` is defined in the `@derived` block using `=` and hence `Free_Ligand` is not used for model fitting but only returned when simulating: │ If `Free_Ligand` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Free_Ligand` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Free_Ligand` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351 ┌ Warning: Variable `Depot_Amount` is defined in the `@derived` block using `=` and hence `Depot_Amount` is not used for model fitting but only returned when simulating: │ If `Depot_Amount` is a random variable, it must be defined in the `@derived` block using `~`; │ if `Depot_Amount` should be returned when simulating, it should be defined in the `@observed` block using `=`; │ if `Depot_Amount` is an intermediate quantity that should not be returned when simulating, it should be defined using `:=`. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/dsl/model_macro.jl:2351
PumasModel
Parameters: tvKa, tvF, tvCl, tvVc, tvKsyn, tvKdeg, tvKSS, tvKint, tvRdepot, tvKSS_depot
Random effects:
Covariates:
Dynamical system variables: Depot, Ltot, Rtot_central
Dynamical system type: Nonlinear ODE
Derived: Free_Ligand, Depot_Amount
Observed: Free_Ligand, Depot_Amount
8.2 SC vs IV Comparison
Show/Hide Code
param_sc = (
tvKa = 0.01,
tvF = 0.6, # SC-specific (ka ~0.24/day, F=60% per Gibiansky)
tvCl = 0.006,
tvVc = 3.0,
tvKsyn = 0.04,
tvKdeg = 0.008,
tvKSS = 1.0,
tvKint = 0.002,
tvRdepot = 2.0,
tvKSS_depot = 1.0, # Depot target binding
)
# SC dosing
pop_sc = [
Subject(
id = 1,
events = DosageRegimen(100, time = 0, cmt = 1), # 100 nmol SC
observations = (Free_Ligand = nothing, Depot_Amount = nothing),
),
]
# IV dosing (use no-FcRn model for comparison)
pop_iv = [
Subject(
id = 1,
events = DosageRegimen(100 * 0.6, time = 0, cmt = 1), # Equivalent absorbed dose (60 nmol)
observations = (Free_Ligand = nothing,),
),
]
sim_sc = simobs(tmdd_sc, pop_sc, param_sc, obstimes = obstimes_fcrn)
sim_iv = simobs(tmdd_no_fcrn, pop_iv, param_no_fcrn, obstimes = obstimes_fcrn)
# Convert to DataFrames
df_iv = DataFrame(sim_iv)
df_iv.Route .= "IV (equivalent dose)"
df_sc = DataFrame(sim_sc)
df_sc.Route .= "SC (with depot binding)"
df_route = vcat(
select(df_iv, :time, :Free_Ligand, :Route),
select(df_sc, :time, :Free_Ligand, :Route),
)
plt =
data(df_route) *
mapping(:time => "Time (hr)", :Free_Ligand => "Free Ligand (nM)"; color = :Route) *
visual(Lines; linewidth = 2)
draw(plt; axis = (yscale = Makie.pseudolog10, title = "IV vs SC Administration"))9 Summary of Special Cases
| Model | Key Feature | When to Use |
|---|---|---|
| Constant \(R_{tot}\) | Receptor pool maintained | Rapid target turnover, large receptor pool |
| Irreversible binding | \(k_{off} = 0\) | Covalent inhibitors, very tight binders |
| Competition | Two drugs, one target | DDI, biosimilar comparison, combination therapy |
| FcRn recycling | Antibody salvage | Therapeutic mAbs, half-life engineering |
| SC with depot binding | Absorption-site target | SC biologics, tissue targets |
Peletier & Gabrielsson (2012) provides rigorous mathematical analysis of TMDD limiting cases, including bounds on receptor concentrations and analytical approximations. Key insights include:
- The minimum free receptor concentration depends on the initial dose and binding parameters
- The time to reach minimum receptor is related to \(1/(k_{on} \cdot L_0)\)
- Recovery kinetics follow receptor synthesis rate \(k_{syn}\)
These analytical results are valuable for understanding special case behavior and validating numerical simulations.
- Start with the simplest model that captures your phenomenon
- Add mechanistic complexity only when scientifically justified
- Consider data availability for parameter estimation
- Validate with independent datasets when possible
For more information on the constant receptor assumption, see Wagner (1976). FcRn biology and its role in mAb pharmacokinetics is reviewed by Roopenian & Akilesh (2007) and Deng et al. (2011). The receptor-mediated endocytosis pathway is discussed by Krippendorff et al. (2009).