TMDD Part 2: The Complete Mechanistic Model

Author

Vijay Ivaturi

1 Introduction

This is Part 2 of our comprehensive TMDD tutorial series. Here we implement the complete mechanistic TMDD model as proposed by Mager & Jusko (2001). This model explicitly tracks ligand, receptor, and complex dynamics.

Tutorial Series Overview:

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

2 Learning Objectives

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

  • Write the complete TMDD differential equations
  • Implement one- and two-compartment full TMDD models in Pumas
  • Understand the role of initial conditions and baseline states
  • Simulate and visualize all model species (ligand, receptor, complex)
  • Handle practical implementation challenges (units, concentrations, stiffness)
  • Perform parameter sensitivity analysis

3 Libraries

using Pumas
using PumasUtilities
using AlgebraOfGraphics
using CairoMakie
using Random
using DataFrames

4 Mathematical Framework

The full TMDD model describes the interaction between a drug (ligand, \(L\)) and its pharmacological target (receptor, \(R\)) through a system of ordinary differential equations.

4.1 One-Compartment Model Equations

For a one-compartment model with IV dosing:

Ligand (free drug) dynamics:

\[ \frac{dL}{dt} = -k_{el} \cdot L - k_{on} \cdot L \cdot R + k_{off} \cdot LR \]

Receptor (free target) dynamics:

\[ \frac{dR}{dt} = k_{syn} - k_{deg} \cdot R - k_{on} \cdot L \cdot R + k_{off} \cdot LR \]

Complex (drug-target) dynamics:

\[ \frac{dLR}{dt} = k_{on} \cdot L \cdot R - k_{off} \cdot LR - k_{int} \cdot LR \]

Where the parameters are:

Parameter Units Description
\(k_{el}\) hr\(^{-1}\) Linear elimination rate constant (\(= CL/V_c\))
\(k_{on}\) nM\(^{-1}\)hr\(^{-1}\) Second-order association rate constant
\(k_{off}\) hr\(^{-1}\) First-order dissociation rate constant
\(k_{syn}\) nM/hr Zero-order receptor synthesis rate
\(k_{deg}\) hr\(^{-1}\) First-order receptor degradation rate
\(k_{int}\) hr\(^{-1}\) First-order complex internalization rate

4.2 Initial Conditions

Before drug administration, the system is at steady state:

  • Ligand: \(L(0) = 0\) (no drug before dosing)
  • Receptor: \(R(0) = R_0 = k_{syn}/k_{deg}\) (baseline receptor level)
  • Complex: \(LR(0) = 0\) (no complex before dosing)

The baseline receptor concentration \(R_0\) is determined by the balance between synthesis and degradation.

4.3 Mass Balance Relationships

The total ligand and total receptor concentrations are conserved:

\[ L_{tot} = L + LR \]

\[ R_{tot} = R + LR \]

These relationships are useful for model verification and data analysis.

5 One-Compartment Full TMDD Model in Pumas

Let’s implement the one-compartment full TMDD model:

tmdd_full_1cmt = @model begin
    @metadata begin
        desc = "One-compartment Full TMDD Model"
        timeu = u"hr"
    end

    @param begin
        """
        Clearance (L/hr)
        """
        tvCl ∈ RealDomain(lower = 0)
        """
        Central volume of distribution (L)
        """
        tvVc ∈ RealDomain(lower = 0)
        """
        Receptor synthesis rate (nM/hr)
        """
        tvKsyn ∈ RealDomain(lower = 0)
        """
        Receptor degradation rate (hr⁻¹)
        """
        tvKdeg ∈ RealDomain(lower = 0)
        """
        Association rate constant (nM⁻¹hr⁻¹)
        """
        tvKon ∈ RealDomain(lower = 0)
        """
        Dissociation rate constant (hr⁻¹)
        """
        tvKoff ∈ RealDomain(lower = 0)
        """
        Internalization rate constant (hr⁻¹)
        """
        tvKint ∈ RealDomain(lower = 0)
    end

    @pre begin
        Cl = tvCl
        Vc = tvVc
        Ksyn = tvKsyn
        Kdeg = tvKdeg
        Kon = tvKon
        Koff = tvKoff
        Kint = tvKint
    end

    @dosecontrol begin
        # Bioavailability trick: dose enters as concentration (amount/Vc)
        bioav = (Ligand = 1 / tvVc,)
    end

    @init begin
        # Baseline receptor level before drug administration
        Receptor = Ksyn / Kdeg
    end

    @vars begin
        # Derived rate constants and quantities
        Kel := Cl / Vc
        R0 := Ksyn / Kdeg
        KD := Koff / Kon

        # Total concentrations (mass balance)
        Ltot := Ligand + Complex
        Rtot := Receptor + Complex

        # Target occupancy and receptor ratio
        TargetOccupancy := Complex / Rtot
        ReceptorRatio := Receptor / R0
    end

    @dynamics begin
        # Free ligand dynamics
        Ligand' = -Kel * Ligand - Kon * Ligand * Receptor + Koff * Complex

        # Free receptor dynamics
        Receptor' = Ksyn - Kdeg * Receptor - Kon * Ligand * Receptor + Koff * Complex

        # Ligand-receptor complex dynamics
        Complex' = Kon * Ligand * Receptor - Koff * Complex - Kint * Complex
    end

    @derived begin
        """
        Free Ligand Concentration (nM)
        """
        Free_Ligand = @. Ligand
        """
        Free Receptor Concentration (nM)
        """
        Free_Receptor = @. Receptor
        """
        Complex Concentration (nM)
        """
        Complex_Conc = @. Complex
        """
        Total Ligand Concentration (nM)
        """
        Total_Ligand = @. Ltot
        """
        Total Receptor Concentration (nM)
        """
        Total_Receptor = @. Rtot
        """
        Target Occupancy (fraction)
        """
        Target_Occupancy = @. TargetOccupancy
    end
end
PumasModel
  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, Target_Occupancy
  Observed: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor, Target_Occupancy
The Bioavailability Trick

Notice the bioav = (Ligand = 1 / tvVc,) in @dosecontrol. This converts the dose (in amount, e.g., mg) directly to concentration (e.g., mg/L or nM) by dividing by the volume of distribution. This approach simplifies the model equations since all state variables are concentrations.

5.1 Setting Up Parameters

# Biologically plausible parameters for a therapeutic mAb
# Based on typical mAb values from the literature
param_1cmt = (
    tvCl = 0.006,     # Clearance: ~0.15 L/day
    tvVc = 3.0,       # Central volume: 3 L (plasma volume)
    tvKsyn = 0.04,    # Synthesis: ~1 nmol/day → R0 = 5 nM
    tvKdeg = 0.008,   # Degradation: ~0.2 day⁻¹ (t½ ~87 hr)
    tvKon = 0.33,     # Association: ~8 nM⁻¹day⁻¹
    tvKoff = 0.33,    # Dissociation: ~8 day⁻¹ (KD = 1 nM)
    tvKint = 0.002,   # Internalization: ~0.04 day⁻¹
)

# Calculate derived quantities
R0 = param_1cmt.tvKsyn / param_1cmt.tvKdeg
KD = param_1cmt.tvKoff / param_1cmt.tvKon
println("Baseline receptor (R0): $(R0) nM")
println("Dissociation constant (KD): $(KD) nM")
Baseline receptor (R0): 5.0 nM
Dissociation constant (KD): 1.0 nM

5.2 Simulating Across Dose Range

Show/Hide Code 
# Create population with multiple doses
# Wide dose range spanning 3 orders of magnitude for clear TMDD visualization
doses = [0.5, 5.0, 50.0, 500.0]  # nmol (equivalent to ~0.075, 0.75, 7.5, 75 mg for a 150 kDa mAb)
dose_labels = ["0.5 nmol", "5 nmol", "50 nmol", "500 nmol"]

Random.seed!(42)

pop_1cmt = [
    Subject(
        id = i,
        events = DosageRegimen(dose, time = 0, cmt = 1),
        covariates = (dose_level = dose,),
        observations = (
            Free_Ligand = nothing,
            Free_Receptor = nothing,
            Complex_Conc = nothing,
            Total_Ligand = nothing,
            Total_Receptor = nothing,
            Target_Occupancy = nothing,
        ),
    ) for (i, dose) in enumerate(doses)
]

# Simulate (extended time for typical mAb kinetics)
obstimes = 0.1:2:672  # Up to 4 weeks

sim_1cmt = simobs(tmdd_full_1cmt, pop_1cmt, param_1cmt, obstimes = obstimes)
Simulated population (Vector{<:Subject})
  Simulated subjects: 4
  Simulated variables: Free_Ligand, Free_Receptor, Complex_Conc, Total_Ligand, Total_Receptor, Target_Occupancy

6 Visualizing Model Output

Show/Hide Code 
# Convert simulation results to DataFrame for plotting
df_1cmt = DataFrame(sim_1cmt)

# Create dose label mapping
dose_label_df = DataFrame(id = string.(1:4), Dose = dose_labels)
df_1cmt = innerjoin(df_1cmt, dose_label_df, on = :id)

6.1 Free Ligand Concentration

Show/Hide Code 
plt =
    data(filter(row -> row.time <= 336, df_1cmt)) *
    mapping(
        :time => "Time (hr)",
        :Free_Ligand => "Free Ligand Concentration (nM)";
        color = :Dose,
    ) *
    visual(Lines; linewidth = 2)

draw(
    plt;
    axis = (yscale = Makie.pseudolog10, title = "One-Compartment Full TMDD: Free Ligand"),
)

Free ligand concentration over time for different doses, showing characteristic TMDD kinetics with dose-dependent terminal slopes.

6.2 Receptor Dynamics

Show/Hide Code 
plt =
    data(filter(row -> row.time <= 504, df_1cmt)) *
    mapping(
        :time => "Time (hr)",
        :Free_Receptor => "Free Receptor Concentration (nM)";
        color = :Dose,
    ) *
    visual(Lines; linewidth = 2)

fg = draw(plt; axis = (title = "One-Compartment Full TMDD: Free Receptor",))

# Add baseline reference line using Makie
hlines!(fg.figure[1, 1], [R0], color = :gray, linestyle = :dash, label = "Baseline (R0)")

fg

Free receptor concentration showing target suppression at higher doses and recovery after drug washout.
Receptor Suppression and Recovery

Higher doses lead to greater suppression of free receptor. After the drug is cleared, receptors recover to baseline through ongoing synthesis. The time course of recovery is governed by \(k_{syn}\) and \(k_{deg}\).

6.3 Complex Formation

Show/Hide Code 
plt =
    data(filter(row -> row.time <= 672, df_1cmt)) *
    mapping(
        :time => "Time (hr)",
        :Complex_Conc => "Complex Concentration (nM)";
        color = :Dose,
    ) *
    visual(Lines; linewidth = 2)

draw(plt; axis = (title = "One-Compartment Full TMDD: Drug-Target Complex",))

Ligand-receptor complex concentration showing dose-dependent formation and eventual clearance via internalization.

6.4 Target Occupancy

Target occupancy is a pharmacodynamically relevant metric showing the fraction of receptors bound by drug:

Show/Hide Code 
plt =
    data(df_1cmt) *
    mapping(
        :time => "Time (hr)",
        :Target_Occupancy => "Target Occupancy (fraction)";
        color = :Dose,
    ) *
    visual(Lines; linewidth = 2)

draw(plt; axis = (title = "One-Compartment Full TMDD: Target Occupancy",))

Target occupancy (fraction of receptors bound) over time. Higher doses maintain occupancy for longer.

6.5 Complete Species Overview

Show/Hide Code 
# Select middle dose for demonstration
dose_idx = 3
df_single = filter(row -> row.id == string(dose_idx), df_1cmt)

# Stack to long format for faceted plot (excluding Target Occupancy which has different y scale)
df_species = stack(
    df_single,
    [:Free_Ligand, :Free_Receptor, :Complex_Conc],
    variable_name = :Species,
    value_name = :Concentration,
)

# Create faceted plot for concentrations
plt =
    data(df_species) *
    mapping(
        :time => "Time (hr)",
        :Concentration => "Concentration (nM)";
        layout = :Species,
    ) *
    visual(Lines; linewidth = 2)

draw(plt; figure = (title = "Dose: $(dose_labels[dose_idx])",))

Complete overview of all species for a single dose (50 nmol), showing the interplay between ligand, receptor, and complex dynamics.

7 Four-Phase Characterization of TMDD Profiles

Peletier & Gabrielsson (2012) provided a detailed characterization of TMDD concentration-time profiles, identifying four distinct phases that occur following drug administration. Understanding these phases is crucial for interpreting PK data and extracting mechanistic information from concentration-time curves.

Deep Dive Available

For a comprehensive exploration of the four phases including parameter identification, constant-rate infusion dynamics, and comparison with the Michaelis-Menten approximation, see Part 1b: TMDD Dynamics Through the Four Phases. That tutorial uses the original Peletier parameters to reproduce key figures from the paper.

7.1 The Four Phases

On a semi-logarithmic plot, the TMDD concentration profile can be divided into four phases:

Phase Name Characteristic Dominant Process
A Initial binding Rapid initial decline Drug-target association
B Target saturation Linear (parallel) decline Linear elimination with saturated target
C Transition Non-linear decline Target-mediated elimination returning
D Terminal Linear (parallel) decline Slowest elimination process

7.2 Phase A: Rapid Initial Binding

Immediately after dosing, free ligand rapidly binds to available receptor. The timescale of this phase is approximately:

\[ \tau_A \sim \frac{1}{k_{on} \cdot L_0} \]

where \(L_0\) is the initial ligand concentration. For high-affinity drugs with large \(k_{on}\), this phase may be very brief and difficult to observe.

7.3 Phase B: Target-Saturated Elimination

Once the target is saturated (\(L \gg K_D\) and \(LR \approx R_{tot}\)), the ligand elimination becomes approximately linear with rate constant:

\[ \lambda_B \approx k_{el} + \frac{k_{int} \cdot R_0}{K_M + L} \]

At high ligand concentrations, this simplifies to approximately \(k_{el}\), representing the non-target-mediated (linear) clearance pathway.

7.4 Phase C: Transition Phase

As ligand concentration decreases and approaches the target concentration, the system transitions from saturated to unsaturated kinetics. This phase exhibits non-linear behavior:

  • Target begins to recover from suppression
  • Target-mediated clearance increases
  • Clearance becomes concentration-dependent

The concentration range for this phase is roughly \(K_D < L < R_0\).

7.5 Phase D: Terminal Phase

At low ligand concentrations (\(L \ll K_D\)), the terminal phase is reached. The terminal slope depends on the relative rates of processes:

\[ \lambda_D \approx k_{int} \text{ (internalization-limited)} \]

or

\[ \lambda_D \approx k_{el} + \frac{k_{int} \cdot R_0}{K_D} \text{ (if linear elimination is slow)} \]

Terminal Slope Interpretation

A key insight from Peletier & Gabrielsson (2012) is that the terminal slope in TMDD is often determined by the internalization rate (\(k_{int}\)) rather than the linear elimination rate (\(k_{el}\)). This differs from the Michaelis-Menten approximation, which can overestimate the terminal elimination rate.

7.6 Visualizing the Four Phases

To clearly visualize all four phases as characterized by Peletier & Gabrielsson (2012), we use parameters that produce distinct phase separation. The key is having:

  • Sufficiently high \(k_{on}\) for visible Phase A
  • Low \(k_{el}\) relative to target-mediated clearance for clear Phase B
  • Appropriate \(R_0\) and \(K_D\) values for visible Phase C transition
  • Slow \(k_{int}\) for extended Phase D
Show/Hide Code 
# Parameters tuned to show clear four-phase behavior
# Based on Peletier & Gabrielsson (2012) characteristic profiles
param_4phase = (
    tvCl = 0.001,      # Low linear clearance to emphasize TMDD
    tvVc = 2.5,        # Central volume
    tvKsyn = 0.5,      # Higher synthesis → R0 = 10 nM
    tvKdeg = 0.05,     # Degradation rate
    tvKon = 1.0,       # Fast association for visible Phase A
    tvKoff = 0.1,      # Slow dissociation (KD = 0.1 nM, high affinity)
    tvKint = 0.02,     # Internalization rate for clear terminal phase
)

R0_4phase = param_4phase.tvKsyn / param_4phase.tvKdeg
KD_4phase = param_4phase.tvKoff / param_4phase.tvKon

# Multiple doses spanning range to show dose-dependent phase behavior
doses_4phase = [10.0, 25.0, 100.0, 250.0, 1000.0]
dose_labels_4phase = ["10 nmol", "25 nmol", "100 nmol", "250 nmol", "1000 nmol"]

pop_4phase = [
    Subject(
        id = i,
        events = DosageRegimen(dose, time = 0, cmt = 1),
        observations = (Free_Ligand = nothing,),
    ) for (i, dose) in enumerate(doses_4phase)
]

# Extended time to capture terminal phase
obstimes_4phase = [0.01; 0.1:0.5:2; 3:1:20; 25:5:100; 110:10:500]

sim_4phase = simobs(tmdd_full_1cmt, pop_4phase, param_4phase, obstimes = obstimes_4phase)
df_4phase = DataFrame(sim_4phase)

# Add dose labels
dose_label_4phase_df = DataFrame(id = string.(1:5), Dose = dose_labels_4phase)
df_4phase = innerjoin(df_4phase, dose_label_4phase_df, on = :id)

# Create the multi-dose plot
fig = Figure(size = (900, 500))
ax = Axis(
    fig[1, 1],
    xlabel = "Time (hr)",
    ylabel = "Free Ligand Concentration (nM)",
    title = "Four-Phase TMDD Profiles Across Doses",
    yscale = log10,
    limits = (nothing, nothing, 0.01, 2000),
)

colors = [:red, :blue, :green, :purple, :orange]
for (i, dose_label) in enumerate(dose_labels_4phase)
    df_dose = filter(row -> row.Dose == dose_label, df_4phase)
    lines!(
        ax,
        df_dose.time,
        max.(df_dose.Free_Ligand, 0.001),
        linewidth = 2,
        color = colors[i],
        label = dose_label,
    )
end

# Add legend
Legend(fig[1, 2], ax, "Dose", framevisible = false)

fig

Four-phase characterization of TMDD profiles across multiple doses, showing the characteristic shape described by Peletier & Gabrielsson (2012).

The plot above demonstrates how the four phases manifest differently across doses:

  • Low doses (10-25 nmol): Phases A and B are brief; the profile quickly transitions through C to D
  • Intermediate doses (100 nmol): All four phases are visible with clear transitions
  • High doses (250-1000 nmol): Extended Phase B with parallel decline; delayed transition to C and D
Show/Hide Code 
# Single dose annotation plot
df_single_4phase = filter(row -> row.Dose == "100 nmol", df_4phase)

fig = Figure(size = (900, 550))
ax = Axis(
    fig[1, 1],
    xlabel = "Time (hr)",
    ylabel = "Free Ligand Concentration (nM)",
    title = "Four-Phase Characterization (100 nmol dose)",
    yscale = log10,
    limits = (nothing, nothing, 0.01, 200),
)

# Plot concentration
lines!(
    ax,
    df_single_4phase.time,
    max.(df_single_4phase.Free_Ligand, 0.001),
    linewidth = 3,
    color = :steelblue,
)

# Phase boundaries (approximate for this parameter set)
vlines!(ax, [1], color = :gray, linestyle = :dash, alpha = 0.6)
vlines!(ax, [15], color = :gray, linestyle = :dash, alpha = 0.6)
vlines!(ax, [80], color = :gray, linestyle = :dash, alpha = 0.6)

# Phase labels with descriptions
text!(ax, 0.3, 80, text = "A", fontsize = 16, font = :bold, color = :darkred)
text!(ax, 5, 25, text = "B", fontsize = 16, font = :bold, color = :darkgreen)
text!(ax, 35, 3, text = "C", fontsize = 16, font = :bold, color = :darkorange)
text!(ax, 200, 0.15, text = "D", fontsize = 16, font = :bold, color = :darkblue)

# Add phase descriptions
text!(ax, 0.15, 25, text = "Rapid\nbinding", fontsize = 10, color = :darkred)
text!(ax, 3, 8, text = "Target\nsaturated", fontsize = 10, color = :darkgreen)
text!(ax, 25, 0.8, text = "Transition", fontsize = 10, color = :darkorange)
text!(ax, 150, 0.05, text = "Terminal", fontsize = 10, color = :darkblue)

# Add reference lines for key concentrations
hlines!(ax, [R0_4phase], color = :gray, linestyle = :dot, alpha = 0.5)
text!(ax, 350, R0_4phase * 1.3, text = "R₀", fontsize = 10, color = :gray)

hlines!(ax, [KD_4phase], color = :gray, linestyle = :dot, alpha = 0.5)
text!(ax, 350, KD_4phase * 1.5, text = "KD", fontsize = 10, color = :gray)

fig

Annotated four-phase profile for a single dose showing phase boundaries and characteristics.

7.7 Parameter Identification from Phases

Each phase provides information about specific TMDD parameters (Peletier & Gabrielsson, 2012):

Phase Observable Parameters Informed
A Duration of rapid decline \(k_{on}\), initial \(L_0\)
B Slope during saturation \(k_{el}\) (linear clearance)
C Curvature, transition point \(K_D\) or \(K_{SS}\), \(R_0\)
D Terminal slope \(k_{int}\) (internalization)
Practical Application

When analyzing TMDD data:

  1. Identify the terminal phase (D) first - its slope gives information about \(k_{int}\)

  2. Look for the transition (C) - the inflection point relates to \(K_D\) and \(R_0\)

  3. The saturated phase (B) slope informs linear clearance \(k_{el}\)

  4. Phase A may require early, dense sampling to observe

    :::

8 Two-Compartment Full TMDD Model

For many drugs, distribution to peripheral tissues is significant. The two-compartment model adds a peripheral compartment for free ligand:

8.1 Equations

The peripheral compartment follows standard linear distribution:

\[ \frac{dL_p}{dt} = k_{12} \cdot L \cdot V_c - k_{21} \cdot L_p \]

And the central ligand equation becomes:

\[ \frac{dL}{dt} = -k_{el} \cdot L - k_{on} \cdot L \cdot R + k_{off} \cdot LR - k_{12} \cdot L + k_{21} \cdot \frac{L_p}{V_c} \]

Where:

  • \(k_{12} = Q/V_c\) (distribution rate from central to peripheral)
  • \(k_{21} = Q/V_p\) (distribution rate from peripheral to central)
  • \(Q\) = Intercompartmental clearance
  • \(V_p\) = Peripheral volume of distribution

8.2 Pumas Implementation

tmdd_full_2cmt = @model begin
    @metadata begin
        desc = "Two-compartment Full TMDD Model"
        timeu = u"hr"
    end

    @param begin
        """
        Clearance (L/hr)
        """
        tvCl ∈ RealDomain(lower = 0)
        """
        Central volume of distribution (L)
        """
        tvVc ∈ RealDomain(lower = 0)
        """
        Intercompartmental clearance (L/hr)
        """
        tvQ ∈ RealDomain(lower = 0)
        """
        Peripheral volume of distribution (L)
        """
        tvVp ∈ RealDomain(lower = 0)
        """
        Receptor synthesis rate (nM/hr)
        """
        tvKsyn ∈ RealDomain(lower = 0)
        """
        Receptor degradation rate (hr⁻¹)
        """
        tvKdeg ∈ RealDomain(lower = 0)
        """
        Association rate constant (nM⁻¹hr⁻¹)
        """
        tvKon ∈ RealDomain(lower = 0)
        """
        Dissociation rate constant (hr⁻¹)
        """
        tvKoff ∈ RealDomain(lower = 0)
        """
        Internalization rate constant (hr⁻¹)
        """
        tvKint ∈ RealDomain(lower = 0)
    end

    @pre begin
        Cl = tvCl
        Vc = tvVc
        Q = tvQ
        Vp = tvVp
        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
        K12 := Q / Vc
        K21 := Q / Vp
        R0 := Ksyn / Kdeg
        KD := Koff / Kon

        Ltot := Ligand + Complex
        Rtot := Receptor + Complex
        TargetOccupancy := Complex / Rtot
    end

    @dynamics begin
        # Free ligand in central compartment
        Ligand' =
            -Kel * Ligand - Kon * Ligand * Receptor + Koff * Complex - K12 * Ligand +
            K21 * (Peripheral / Vc)

        # Free ligand in peripheral compartment (as amount)
        Peripheral' = K12 * Ligand * Vc - K21 * Peripheral

        # Free receptor dynamics (central only)
        Receptor' = Ksyn - Kdeg * Receptor - Kon * Ligand * Receptor + Koff * Complex

        # Complex dynamics (central only)
        Complex' = Kon * Ligand * Receptor - Koff * Complex - Kint * Complex
    end

    @derived begin
        Free_Ligand = @. Ligand
        Free_Receptor = @. Receptor
        Complex_Conc = @. Complex
        Peripheral_Conc = @. Peripheral / Vp
        Total_Ligand = @. Ltot
        Target_Occupancy = @. TargetOccupancy
    end
end
PumasModel
  Parameters: tvCl, tvVc, tvQ, tvVp, tvKsyn, tvKdeg, tvKon, tvKoff, tvKint
  Random effects:
  Covariates:
  Dynamical system variables: Ligand, Peripheral, Receptor, Complex
  Dynamical system type: Nonlinear ODE
  Derived: Free_Ligand, Free_Receptor, Complex_Conc, Peripheral_Conc, Total_Ligand, Target_Occupancy
  Observed: Free_Ligand, Free_Receptor, Complex_Conc, Peripheral_Conc, Total_Ligand, Target_Occupancy

8.3 Simulating the Two-Compartment Model

Show/Hide Code 
# Parameters for two-compartment model
# Based on typical mAb values from the literature
param_2cmt = (
    tvCl = 0.006,     # Clearance: ~0.15 L/day
    tvVc = 3.0,       # Central volume: 3 L
    tvQ = 0.02,       # Intercompartmental clearance: ~0.45 L/day
    tvVp = 3.0,       # Peripheral volume: 3 L
    tvKsyn = 0.04,    # Synthesis rate: ~1 nmol/day
    tvKdeg = 0.008,   # Degradation rate: ~0.2 day⁻¹
    tvKon = 0.33,     # Association rate: ~8 nM⁻¹day⁻¹
    tvKoff = 0.33,    # Dissociation rate: ~8 day⁻¹
    tvKint = 0.002,   # Internalization rate: ~0.04 day⁻¹
)

pop_2cmt = [
    Subject(
        id = i,
        events = DosageRegimen(dose, time = 0, cmt = 1),
        observations = (
            Free_Ligand = nothing,
            Peripheral_Conc = nothing,
            Free_Receptor = nothing,
            Complex_Conc = nothing,
        ),
    ) for (i, dose) in enumerate(doses)
]

sim_2cmt = simobs(tmdd_full_2cmt, pop_2cmt, param_2cmt, obstimes = obstimes)
Simulated population (Vector{<:Subject})
  Simulated subjects: 4
  Simulated variables: Free_Ligand, Free_Receptor, Complex_Conc, Peripheral_Conc, Total_Ligand, Target_Occupancy

8.4 Comparing One- and Two-Compartment Models

Show/Hide Code 
dose_idx = 3  # 10 mg/kg

# Prepare data for model comparison
df_2cmt = DataFrame(sim_2cmt)
df_2cmt = innerjoin(df_2cmt, dose_label_df, on = :id)

# Filter to dose of interest and early time
df_1cmt_sub = filter(row -> row.id == string(dose_idx) && row.time <= 24, df_1cmt)
df_2cmt_sub = filter(row -> row.id == string(dose_idx) && row.time <= 24, df_2cmt)

# Panel 1: Model comparison
df_1cmt_sub.Model .= "1-compartment"
df_2cmt_sub.Model .= "2-compartment"
df_comp = vcat(
    select(df_1cmt_sub, :time, :Free_Ligand, :Model),
    select(df_2cmt_sub, :time, :Free_Ligand, :Model),
)

plt1 =
    data(df_comp) *
    mapping(
        :time => "Time (hr)",
        :Free_Ligand => "Free Ligand (nM)";
        color = :Model,
        linestyle = :Model,
    ) *
    visual(Lines; linewidth = 2)

# Panel 2: Central vs Peripheral for 2-cmt
df_cmt_long = stack(
    df_2cmt_sub,
    [:Free_Ligand, :Peripheral_Conc],
    variable_name = :Compartment,
    value_name = :Concentration,
)
df_cmt_long.Compartment =
    replace.(
        string.(df_cmt_long.Compartment),
        "Free_Ligand" => "Central",
        "Peripheral_Conc" => "Peripheral",
    )

plt2 =
    data(df_cmt_long) *
    mapping(
        :time => "Time (hr)",
        :Concentration => "Concentration (nM)";
        color = :Compartment,
        linestyle = :Compartment,
    ) *
    visual(Lines; linewidth = 2)

fig = Figure(size = (900, 400))
draw!(
    fig[1, 1],
    plt1;
    axis = (yscale = Makie.pseudolog10, title = "Free Ligand Comparison"),
)
draw!(
    fig[1, 2],
    plt2;
    axis = (yscale = Makie.pseudolog10, title = "Central vs Peripheral (2-cmt)"),
)
fig

Comparison of one- and two-compartment full TMDD models. The two-compartment model shows initial distribution phase.

9 Parameter Sensitivity Analysis

Understanding how each parameter affects model output is crucial for model development and interpretation.

Show/Hide Code 
# Base parameter set (literature-based values)
base_params = (
    tvCl = 0.006,
    tvVc = 3.0,
    tvKsyn = 0.04,
    tvKdeg = 0.008,
    tvKon = 0.33,
    tvKoff = 0.33,
    tvKint = 0.002,
)

# Single subject at fixed dose
single_pop = [
    Subject(
        id = 1,
        events = DosageRegimen(50, time = 0, cmt = 1),
        observations = (Free_Ligand = nothing,),
    ),
]

# Parameters to vary and their ranges
param_variations = [
    (:tvKon, [0.1, 0.33, 0.5, 1.0], "kon"),
    (:tvKoff, [0.1, 0.33, 1.0, 3.0], "koff"),
    (:tvKint, [0.001, 0.002, 0.005, 0.01], "kint"),
    (:tvKdeg, [0.004, 0.008, 0.016, 0.032], "kdeg"),
]

# Build DataFrame with all parameter variations
df_sens = DataFrame()
for (param_name, values, label) in param_variations
    for val in values
        modified_params = merge(base_params, NamedTuple{(param_name,)}((val,)))
        sim = simobs(tmdd_full_1cmt, single_pop, modified_params, obstimes = obstimes)
        sim_df = DataFrame(sim)
        sim_df.Parameter .= label
        sim_df.Value .= string(val)
        append!(df_sens, sim_df)
    end
end

plt =
    data(df_sens) *
    mapping(
        :time => "Time (hr)",
        :Free_Ligand => "Free Ligand (nM)";
        color = :Value => "Value",
        layout = :Parameter,
    ) *
    visual(Lines; linewidth = 2)

draw(plt; axis = (yscale = Makie.pseudolog10,))

Sensitivity of free ligand concentration to key TMDD parameters. Each panel varies one parameter while holding others constant.

10 Practical Implementation Considerations

10.1 Units and Scaling

When implementing TMDD models, careful attention to units is essential:

Quantity Common Units Alternative
Concentration nM, pM ng/mL, μg/mL
Volume L mL
Clearance L/hr mL/min
Rate constants hr⁻¹ day⁻¹, min⁻¹
\(k_{on}\) nM⁻¹hr⁻¹ M⁻¹s⁻¹
Unit Consistency

Ensure all parameters have consistent units. A common source of error is mixing concentration units (e.g., nM for ligand, ng/mL for parameters).

10.2 Handling Stiff ODEs

The full TMDD model can be numerically stiff, especially when:

  • \(k_{on}\) is large (fast binding)
  • There are large concentration gradients

Pumas uses adaptive ODE solvers that handle stiffness automatically. For particularly challenging systems, you can specify solver options.

10.3 Initial Condition Considerations

The receptor initial condition \(R_0 = k_{syn}/k_{deg}\) assumes steady state before dosing. If this assumption doesn’t hold (e.g., disease state alters baseline), you can:

  1. Estimate \(R_0\) directly as a parameter
  2. Model disease-related changes in \(k_{syn}\) or \(k_{deg}\)
  3. Use time-varying parameters
# Example: Estimating R0 directly instead of deriving from ksyn/kdeg
tmdd_with_R0 = @model begin
    @param begin
        tvCl ∈ RealDomain(lower = 0)
        tvVc ∈ RealDomain(lower = 0)
        tvR0 ∈ RealDomain(lower = 0)      # Directly estimated
        tvKdeg ∈ RealDomain(lower = 0)
        tvKon ∈ RealDomain(lower = 0)
        tvKoff ∈ RealDomain(lower = 0)
        tvKint ∈ RealDomain(lower = 0)
    end

    @pre begin
        Cl = tvCl
        Vc = tvVc
        R0 = tvR0
        Kdeg = tvKdeg
        Kon = tvKon
        Koff = tvKoff
        Kint = tvKint
        # Derive Ksyn to maintain consistency
        Ksyn = R0 * Kdeg
    end

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

    @init begin
        Receptor = R0
    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
    end
end
PumasModel
  Parameters: tvCl, tvVc, tvR0, tvKdeg, tvKon, tvKoff, tvKint
  Random effects:
  Covariates:
  Dynamical system variables: Ligand, Receptor, Complex
  Dynamical system type: Nonlinear ODE
  Derived: Free_Ligand
  Observed: Free_Ligand

11 Verifying Model Behavior

11.1 Mass Balance Check

A useful verification is to check that mass balance holds:

Show/Hide Code 
# Verify mass balance: Ltot should be conserved (minus elimination)
df_mass = filter(row -> row.id == "3", df_1cmt)

plt =
    data(df_mass) *
    mapping(:time => "Time (hr)", :Total_Ligand => "Total Ligand (nM)") *
    visual(Lines; linewidth = 2)

draw(plt; axis = (title = "Mass Balance Verification: Total Ligand Over Time",))

11.2 Steady-State Verification

Before dosing, the receptor should be at steady state:

Show/Hide Code 
# Simulate without dosing to verify steady state
no_dose_pop = [Subject(id = 1, observations = (Free_Receptor = nothing,))]

sim_nodose = simobs(tmdd_full_1cmt, no_dose_pop, param_1cmt, obstimes = 0:10:200)

println("Expected R0: $(R0) nM")
println("Simulated R (should be constant at R0):")
println("  t=0: $(sim_nodose[1].observations.Free_Receptor[1]) nM")
println("  t=100: $(sim_nodose[1].observations.Free_Receptor[11]) nM")
println("  t=200: $(sim_nodose[1].observations.Free_Receptor[21]) nM")
Expected R0: 5.0 nM
Simulated R (should be constant at R0):
  t=0: 5.0 nM
  t=100: 5.0 nM
  t=200: 5.0 nM

12 Summary

In this tutorial, we covered:

  • Mathematical framework: Complete ODE system for full TMDD
  • One-compartment model: Basic implementation with ligand, receptor, and complex
  • Two-compartment model: Extension with peripheral distribution
  • Visualization: All species dynamics across dose range
  • Parameter sensitivity: How each parameter affects model output
  • Practical considerations: Units, stiffness, initial conditions

The full TMDD model provides the most mechanistic description but may be overparameterized for typical data. In Part 3, we explore simplified approximations (QE, QSS, MM) that maintain key TMDD features while improving identifiability.

The full TMDD model was first introduced by Mager & Jusko (2001), with subsequent refinements by Mager & Krzyzanski (2005). The four-phase characterization of TMDD profiles was described in detail by Peletier & Gabrielsson (2012), and comprehensive tutorials are available from Dua et al. (2015).

13 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
Mager, D. E., & Jusko, W. J. (2001). General pharmacokinetic model for drugs exhibiting target-mediated drug disposition. Journal of Pharmacokinetics and Pharmacodynamics, 28(6), 507–532. https://doi.org/10.1023/A:1014414520282
Mager, D. E., & Krzyzanski, W. (2005). Quasi-equilibrium pharmacokinetic model for drugs exhibiting target-mediated drug disposition. Pharmaceutical Research, 22(10), 1589–1596. https://doi.org/10.1007/s11095-005-6650-0
Peletier, L. A., & Gabrielsson, J. (2012). Dynamics of target-mediated drug disposition: Characteristic profiles and parameter identification. Journal of Pharmacokinetics and Pharmacodynamics, 39(5), 429–451. https://doi.org/10.1007/s10928-012-9260-6

Reuse