Indirect Response Pharmacodynamic Models

Authors

Vijay Ivaturi

Marwa Elsaeed

1 Learning Objectives

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

Understand the turnover (indirect response) framework
Distinguish between the four indirect response model types
Choose the appropriate IDR type based on drug mechanism
Implement all four IDR models in Pumas
Compare response characteristics across model types
Apply IDR models to real-world scenarios

2 Introduction

Indirect response (IDR) models describe drug effects that are mediated through modulation of physiological processes with inherent production and elimination rates (Dayneka, Garg, and Jusko 1993; Sharma and Jusko 1998). Unlike direct response models, where effect immediately follows concentration, IDR models capture situations where:

Effect persists after drug concentration declines
Drug action is mediated through physiological turnover
Response involves synthesis/degradation of endogenous substances

flowchart TD
    subgraph Physiological System
        kin[k_in<br/>Production Rate] --> R[Response<br/>Variable R]
        R --> kout[k_out<br/>Removal Rate]
    end

    Drug --> |Inhibit or<br/>Stimulate| kin
    Drug --> |Inhibit or<br/>Stimulate| kout

    style R fill:#e8f5e9

Figure 1: Indirect response model: Drug affects the balance between production and removal

2.1 The Turnover Concept

Many physiological variables are maintained at homeostasis through a balance between production (k_{in}) and removal (k_{out}):

\frac{dR}{dt} = k_{in} - k_{out} \cdot R

At steady state (baseline, no drug):

\frac{dR}{dt} = 0 \implies R_{ss} = \frac{k_{in}}{k_{out}}

Key Insight

The drug doesn’t directly cause the effect. Instead, it modifies the rate of production or removal, and the effect emerges as the system moves to a new equilibrium. This is why effects often outlast drug concentrations.

2.2 Clinical Examples of Turnover Processes

Response Variable	Production	Removal	Turnover Time
Clotting factors	Hepatic synthesis	Degradation	1-3 days
Blood glucose	Hepatic production	Cellular uptake	Hours
Cortisol	Adrenal secretion	Metabolism	~1 hour
Platelets	Bone marrow	Clearance	7-10 days
LDL cholesterol	Hepatic synthesis	LDL receptors	~3 days

3 The Four Indirect Response Types

Drugs can affect turnover processes in four distinct ways:

flowchart LR
    subgraph Type1["IDR Type I"]
        direction TB
        kin1[k_in] -->|Drug INHIBITS| R1[R decreases]
    end

    subgraph Type2["IDR Type II"]
        direction TB
        kout2[k_out] -->|Drug INHIBITS| R2[R increases]
    end

    subgraph Type3["IDR Type III"]
        direction TB
        kin3[k_in] -->|Drug STIMULATES| R3[R increases]
    end

    subgraph Type4["IDR Type IV"]
        direction TB
        kout4[k_out] -->|Drug STIMULATES| R4[R decreases]
    end

    style Type1 fill:#ffebee
    style Type2 fill:#e8f5e9
    style Type3 fill:#e8f5e9
    style Type4 fill:#ffebee

Figure 2: The four types of indirect response models

Type	Drug Action	Effect on k_{in}	Effect on k_{out}	Response Direction
I	Inhibition of k_{in}	↓	—	Decreases
II	Inhibition of k_{out}	—	↓	Increases
III	Stimulation of k_{in}	↑	—	Increases
IV	Stimulation of k_{out}	—	↑	Decreases

4 Mathematical Framework

4.1 General IDR Equations

Baseline condition (no drug):

\frac{dR}{dt} = k_{in,0} - k_{out} \cdot R

With drug effect, the modulated rate becomes:

\frac{dR}{dt} = k_{in,0} \cdot \left(1 - \frac{I_{max} \cdot C}{IC_{50} + C}\right) - k_{out} \cdot R

The production rate is reduced by the drug.

\frac{dR}{dt} = k_{in,0} - k_{out} \cdot \left(1 - \frac{I_{max} \cdot C}{IC_{50} + C}\right) \cdot R

The removal rate is reduced by the drug.

\frac{dR}{dt} = k_{in,0} \cdot \left(1 + \frac{E_{max} \cdot C}{EC_{50} + C}\right) - k_{out} \cdot R

The production rate is increased by the drug.

\frac{dR}{dt} = k_{in,0} - k_{out} \cdot \left(1 + \frac{E_{max} \cdot C}{EC_{50} + C}\right) \cdot R

The removal rate is increased by the drug.

4.2 Key Parameters

Parameter	Description	Units
k_{in,0}	Baseline production rate	response/time
k_{out}	First-order removal rate constant	1/time
I_{max}	Maximum fractional inhibition (0 to 1)	dimensionless
E_{max}	Maximum fractional stimulation	dimensionless
IC_{50}/EC_{50}	Concentration for 50% effect	mass/volume

Turnover Half-Life

The turnover half-life determines how quickly the response variable returns to baseline after drug removal:

t_{1/2,turnover} = \frac{\ln(2)}{k_{out}}

A process with t_{1/2} = 24 hours (e.g., clotting factors) takes ~4-5 days to fully recover after drug cessation.

5 Pumas Implementation

5.1 IDR Type I: Inhibition of Production

idr_type1 = @model begin
    @metadata begin
        desc = "Indirect Response Model Type I: Inhibition of kin"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
        pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
        pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)

        # PD parameters
        """
        Baseline response (R0 = kin0/kout)
        """
        pop_R0 ∈ RealDomain(; lower = 0, init = 100.0)
        """
        Turnover time (1/kout)
        """
        pop_turnover ∈ RealDomain(; lower = 0, init = 10.0)
        """
        Maximum fractional inhibition of kin
        """
        pop_Imax ∈ RealDomain(; lower = 0, upper = 1, init = 0.9)
        """
        Concentration at 50% inhibition
        """
        pop_IC50 ∈ RealDomain(; lower = 0, init = 2.0)

        # BSV and RUV
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(1)
        σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add_pd ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        # PK parameters
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])

        # PD parameters
        R0 = pop_R0 * exp(η_pd[1])
        kout = 1 / pop_turnover
        kin0 = R0 * kout  # Derived from steady-state condition
        Imax = pop_Imax
        IC50 = pop_IC50
    end

    @init begin
        Response = R0  # Start at baseline
    end

    @vars begin
        # Plasma concentration
        Cp = Central / Vc
        # IDR Type I: Drug inhibits production
        inhibition = Imax * Cp / (IC50 + Cp)
        kin = kin0 * (1 - inhibition)
        # Response for plotting
        Resp = Response
    end

    @dynamics begin
        # PK
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
        # Response dynamics
        Response' = kin - kout * Response
    end

    @derived begin
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        response ~ @. Normal(Response, σ_add_pd)
    end
end

PumasModel
  Parameters: pop_Ka, pop_CL, pop_Vc, pop_R0, pop_turnover, pop_Imax, pop_IC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
  Random effects: η_pk, η_pd
  Covariates:
  Dynamical system variables: Depot, Central, Response
  Dynamical system type: Nonlinear ODE
  Derived: conc, response, Cp, inhibition, kin, Resp
  Observed: conc, response, Cp, inhibition, kin, Resp

5.2 IDR Type II: Inhibition of Removal

idr_type2 = @model begin
    @metadata begin
        desc = "Indirect Response Model Type II: Inhibition of kout"
        timeu = u"hr"
    end

    @param begin
        pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
        pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
        pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)
        pop_R0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_turnover ∈ RealDomain(; lower = 0, init = 10.0)
        pop_Imax ∈ RealDomain(; lower = 0, upper = 1, init = 0.9)
        pop_IC50 ∈ RealDomain(; lower = 0, init = 2.0)
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(1)
        σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add_pd ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])
        R0 = pop_R0 * exp(η_pd[1])
        kout0 = 1 / pop_turnover
        kin = R0 * kout0
        Imax = pop_Imax
        IC50 = pop_IC50
    end

    @init begin
        Response = R0
    end

    @vars begin
        # Plasma concentration
        Cp = Central / Vc
        # IDR Type II: Drug inhibits removal
        inhibition = Imax * Cp / (IC50 + Cp)
        kout = kout0 * (1 - inhibition)
        # Response for plotting
        Resp = Response
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
        # Response dynamics
        Response' = kin - kout * Response
    end

    @derived begin
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        response ~ @. Normal(Response, σ_add_pd)
    end
end

PumasModel
  Parameters: pop_Ka, pop_CL, pop_Vc, pop_R0, pop_turnover, pop_Imax, pop_IC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
  Random effects: η_pk, η_pd
  Covariates:
  Dynamical system variables: Depot, Central, Response
  Dynamical system type: Nonlinear ODE
  Derived: conc, response, Cp, inhibition, kout, Resp
  Observed: conc, response, Cp, inhibition, kout, Resp

5.3 IDR Type III: Stimulation of Production

idr_type3 = @model begin
    @metadata begin
        desc = "Indirect Response Model Type III: Stimulation of kin"
        timeu = u"hr"
    end

    @param begin
        pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
        pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
        pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)
        pop_R0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_turnover ∈ RealDomain(; lower = 0, init = 10.0)
        """
        Maximum fractional stimulation of kin
        """
        pop_Emax ∈ RealDomain(; lower = 0, init = 2.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 2.0)
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(1)
        σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add_pd ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])
        R0 = pop_R0 * exp(η_pd[1])
        kout = 1 / pop_turnover
        kin0 = R0 * kout
        Emax = pop_Emax
        EC50 = pop_EC50
    end

    @init begin
        Response = R0
    end

    @vars begin
        # Plasma concentration
        Cp = Central / Vc
        # IDR Type III: Drug stimulates production
        stimulation = Emax * Cp / (EC50 + Cp)
        kin = kin0 * (1 + stimulation)
        # Response for plotting
        Resp = Response
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
        # Response dynamics
        Response' = kin - kout * Response
    end

    @derived begin
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        response ~ @. Normal(Response, σ_add_pd)
    end
end

PumasModel
  Parameters: pop_Ka, pop_CL, pop_Vc, pop_R0, pop_turnover, pop_Emax, pop_EC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
  Random effects: η_pk, η_pd
  Covariates:
  Dynamical system variables: Depot, Central, Response
  Dynamical system type: Nonlinear ODE
  Derived: conc, response, Cp, stimulation, kin, Resp
  Observed: conc, response, Cp, stimulation, kin, Resp

5.4 IDR Type IV: Stimulation of Removal

idr_type4 = @model begin
    @metadata begin
        desc = "Indirect Response Model Type IV: Stimulation of kout"
        timeu = u"hr"
    end

    @param begin
        pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
        pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
        pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)
        pop_R0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_turnover ∈ RealDomain(; lower = 0, init = 10.0)
        pop_Emax ∈ RealDomain(; lower = 0, init = 2.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 2.0)
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(1)
        σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add_pd ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])
        R0 = pop_R0 * exp(η_pd[1])
        kout0 = 1 / pop_turnover
        kin = R0 * kout0
        Emax = pop_Emax
        EC50 = pop_EC50
    end

    @init begin
        Response = R0
    end

    @vars begin
        # Plasma concentration
        Cp = Central / Vc
        # IDR Type IV: Drug stimulates removal
        stimulation = Emax * Cp / (EC50 + Cp)
        kout = kout0 * (1 + stimulation)
        # Response for plotting
        Resp = Response
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
        # Response dynamics
        Response' = kin - kout * Response
    end

    @derived begin
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        response ~ @. Normal(Response, σ_add_pd)
    end
end

PumasModel
  Parameters: pop_Ka, pop_CL, pop_Vc, pop_R0, pop_turnover, pop_Emax, pop_EC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
  Random effects: η_pk, η_pd
  Covariates:
  Dynamical system variables: Depot, Central, Response
  Dynamical system type: Nonlinear ODE
  Derived: conc, response, Cp, stimulation, kout, Resp
  Observed: conc, response, Cp, stimulation, kout, Resp

6 Comparative Simulation

Let’s simulate all four IDR types with the same PK parameters to compare their behavior:

# Parameters for inhibition models (Types I, II)
params_inhib = (
    pop_Ka = 1.0,
    pop_CL = 0.5,
    pop_Vc = 20.0,
    pop_R0 = 100.0,
    pop_turnover = 10.0,  # t1/2 ≈ 7 hr
    pop_Imax = 0.8,
    pop_IC50 = 2.0,
    Ω_pk = Diagonal([0.0, 0.0, 0.0]),  # No variability for comparison
    Ω_pd = Diagonal([0.0]),
    σ_prop_pk = 0.0,
    σ_add_pd = 0.0,
)

# Parameters for stimulation models (Types III, IV)
params_stim = (
    pop_Ka = 1.0,
    pop_CL = 0.5,
    pop_Vc = 20.0,
    pop_R0 = 100.0,
    pop_turnover = 10.0,  # t1/2 ≈ 7 hr
    pop_Emax = 1.5,
    pop_EC50 = 2.0,
    Ω_pk = Diagonal([0.0, 0.0, 0.0]),  # No variability for comparison
    Ω_pd = Diagonal([0.0]),
    σ_prop_pk = 0.0,
    σ_add_pd = 0.0,
)

# Dosing and subject
dosing = DosageRegimen(100, time = 0, cmt = 1)
subject = Subject(id = 1, events = dosing)

# Simulate all four types (no residual error for clean visualization)
sim1 = simobs(idr_type1, subject, params_inhib, obstimes = 0:0.5:72, simulate_error = false)
sim2 = simobs(idr_type2, subject, params_inhib, obstimes = 0:0.5:72, simulate_error = false)
sim3 = simobs(idr_type3, subject, params_stim, obstimes = 0:0.5:72, simulate_error = false)
sim4 = simobs(idr_type4, subject, params_stim, obstimes = 0:0.5:72, simulate_error = false)

SimulatedObservations
  Simulated variables: conc, response, Cp, stimulation, kout, Resp
  Time: 0.0:0.5:72.0

# Convert simulations to DataFrames
df1 = DataFrame(sim1)
df1.idr_type .= "Type I: Inhibition of kin"
df2 = DataFrame(sim2)
df2.idr_type .= "Type II: Inhibition of kout"
df3 = DataFrame(sim3)
df3.idr_type .= "Type III: Stimulation of kin"
df4 = DataFrame(sim4)
df4.idr_type .= "Type IV: Stimulation of kout"

# PK plot (use df1 since PK is same for all) - using vars: Cp
pk_plt =
    data(df1) *
    mapping(:time => "Time (hr)", :Cp => "Concentration (mg/L)") *
    visual(Lines, linewidth = 2, color = :steelblue)

fig = Figure(size = (900, 700))
draw!(fig[1, 1:2], pk_plt; axis = (title = "Plasma Concentration (Common to All)",))

# Draw individual response panels
for (i, (df, title, color, pos)) in enumerate(
    zip(
        [df1, df2, df3, df4],
        [
            "Type I: Inhibition of kin",
            "Type II: Inhibition of kout",
            "Type III: Stimulation of kin",
            "Type IV: Stimulation of kout",
        ],
        [:red, :green, :purple, :orange],
        [(2, 1), (2, 2), (3, 1), (3, 2)],
    ),
)
    plt =
        data(df) *
        mapping(:time => "Time (hr)", :Resp => "Response") *
        visual(Lines, linewidth = 2, color = color)
    ax = Axis(fig[pos...], xlabel = "Time (hr)", ylabel = "Response", title = title)
    draw!(ax, plt)
    hlines!(ax, [100], linestyle = :dash, color = :gray)
end

fig

Figure 3: Comparison of all four IDR model types with identical PK

6.1 Key Observations

Response Characteristics

Types I and IV → Response decreases

Type I: Reduced production leads to gradual decline
Type IV: Increased removal leads to faster decline

Types II and III → Response increases

Type II: Reduced removal causes accumulation
Type III: Increased production causes rise

All types: Response returns to baseline after drug clears, governed by the turnover half-life.

6.2 Overlay Comparison

# Prepare labeled DataFrames (reuse df1-df4 from previous block) - using vars: Resp
df1_labeled = select(df1, :time, :Resp)
df1_labeled.type .= "Type I (↓kin)"
df2_labeled = select(df2, :time, :Resp)
df2_labeled.type .= "Type II (↓kout)"
df3_labeled = select(df3, :time, :Resp)
df3_labeled.type .= "Type III (↑kin)"
df4_labeled = select(df4, :time, :Resp)
df4_labeled.type .= "Type IV (↑kout)"

overlay_df = vcat(df1_labeled, df2_labeled, df3_labeled, df4_labeled)

plt =
    data(overlay_df) *
    mapping(:time => "Time (hr)", :Resp => "Response", color = :type => "IDR Type") *
    visual(Lines, linewidth = 2)

draw(
    plt;
    axis = (
        xlabel = "Time (hr)",
        ylabel = "Response",
        title = "All Four IDR Types Compared",
        width = 700,
        height = 400,
    ),
)

Figure 4: Overlaid comparison of all four IDR types

7 Choosing the Right IDR Type

The choice of IDR model depends on the mechanism of drug action:

Inhibition of production rate

Drug	Target	Response Variable
Corticosteroids	ACTH release	Cortisol levels
Statins	HMG-CoA reductase	LDL cholesterol
Methotrexate	DNA synthesis	Cell proliferation

Clinical scenario: A drug suppresses the synthesis of a biomarker. After drug removal, the biomarker slowly returns as production resumes (Chakraborty, Krzyzanski, and Jusko 1999).

Inhibition of removal rate

Drug	Target	Response Variable
PPIs	H⁺/K⁺-ATPase	Gastric acid
ACE inhibitors	Angiotensin degradation	Bradykinin levels
Warfarin	Clotting factor degradation	Clotting factors

Clinical scenario: A drug blocks the clearance of a substance, causing it to accumulate above baseline (Hamberg et al. 2007).

Stimulation of production rate

Drug	Target	Response Variable
Erythropoietin	RBC production	Hemoglobin
G-CSF	Neutrophil production	WBC count
Insulin secretagogues	β-cell secretion	Insulin levels

Clinical scenario: A drug enhances the synthesis of an endogenous factor, leading to supranormal levels during treatment.

Stimulation of removal rate

Drug	Target	Response Variable
Diuretics	Sodium excretion	Body sodium
Enzyme inducers	Metabolic clearance	Drug levels
Laxatives	GI motility	Gut contents

Clinical scenario: A drug accelerates the elimination of a substance, causing levels to drop below baseline.

8 Effect of Turnover Time

The turnover time (1/k_{out}) significantly affects the response dynamics:

turnover_times = [2.0, 10.0, 30.0]
labels = ["τ = 2 hr (fast)", "τ = 10 hr (moderate)", "τ = 30 hr (slow)"]

# Simulate for each turnover time and collect results (no residual error for visualization)
turnover_results = map(zip(turnover_times, labels)) do (τ, label)
    params_var = merge(params_inhib, (; pop_turnover = τ))
    sim_var = simobs(
        idr_type1,
        subject,
        params_var,
        obstimes = 0:0.5:120,
        simulate_error = false,
    )
    df = DataFrame(sim_var)
    df.turnover_label .= label
    df
end
turnover_df = vcat(turnover_results...)

plt =
    data(turnover_df) *
    mapping(
        :time => "Time (hr)",
        :Resp => "Response",
        color = :turnover_label => "Turnover Time",
    ) *
    visual(Lines, linewidth = 2)

draw(
    plt;
    axis = (
        xlabel = "Time (hr)",
        ylabel = "Response",
        title = "Impact of Turnover Time on Response",
        width = 700,
        height = 350,
    ),
)

Figure 5: Effect of turnover time on response dynamics (Type I model)

Interpretation

Fast turnover (τ = 2 hr): Response changes quickly, returns to baseline fast
Slow turnover (τ = 30 hr): Response changes gradually, takes days to normalize

The turnover time determines the “memory” of the system.

9 Fitting Workflow (Type I Example)

9.1 Generate Synthetic Data

# Parameters with variability
params_fit = (
    pop_Ka = 1.0,
    pop_CL = 0.5,
    pop_Vc = 20.0,
    pop_R0 = 100.0,
    pop_turnover = 10.0,
    pop_Imax = 0.8,
    pop_IC50 = 2.0,
    Ω_pk = Diagonal([0.09, 0.09, 0.09]),
    Ω_pd = Diagonal([0.04]),
    σ_prop_pk = 0.1,
    σ_add_pd = 5.0,
)

# Population
population = map(1:20) do i
    Subject(id = i, events = dosing)
end

Random.seed!(789)
sim_data =
    simobs(idr_type1, population, params_fit, obstimes = [1, 2, 4, 8, 12, 24, 36, 48])
df_sim = DataFrame(sim_data)

9.2 Fit the Model

fit_pop = read_pumas(
    df_sim;
    id = :id,
    time = :time,
    amt = :amt,
    evid = :evid,
    cmt = :cmt,
    observations = [:conc, :response],
)

init_params = (
    pop_Ka = 0.8,
    pop_CL = 0.4,
    pop_Vc = 25.0,
    pop_R0 = 90.0,
    pop_turnover = 8.0,
    pop_Imax = 0.7,
    pop_IC50 = 1.5,
    Ω_pk = Diagonal([0.05, 0.05, 0.05]),
    Ω_pd = Diagonal([0.02]),
    σ_prop_pk = 0.15,
    σ_add_pd = 8.0,
)

fit_result = fit(idr_type1, fit_pop, init_params, FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     7.043609e+02     1.362453e+02
 * time: 0.04104495048522949
     1     6.660506e+02     9.658899e+01
 * time: 2.253674030303955
     2     6.381392e+02     2.852823e+01
 * time: 2.664001941680908
     3     6.325945e+02     1.932805e+01
 * time: 3.0173418521881104
     4     6.318110e+02     2.025626e+01
 * time: 3.312253952026367
     5     6.305202e+02     1.783273e+01
 * time: 3.598180055618286
     6     6.296707e+02     2.549335e+01
 * time: 3.8941099643707275
     7     6.276822e+02     1.160172e+01
 * time: 4.17707896232605
     8     6.262329e+02     6.178991e+00
 * time: 4.470814943313599
     9     6.259682e+02     4.909688e+00
 * time: 4.751853942871094
    10     6.258816e+02     1.990152e+00
 * time: 5.04046106338501
    11     6.258309e+02     1.195364e+00
 * time: 5.31663703918457
    12     6.257060e+02     1.646875e+00
 * time: 5.604434967041016
    13     6.254642e+02     1.535593e+00
 * time: 5.887721061706543
    14     6.252608e+02     2.334304e+00
 * time: 6.178802967071533
    15     6.252488e+02     3.811302e+00
 * time: 6.458811044692993
    16     6.252354e+02     1.398465e+00
 * time: 6.741858005523682
    17     6.252320e+02     3.406503e-01
 * time: 7.017269849777222
    18     6.252302e+02     3.004677e-01
 * time: 7.29183292388916
    19     6.252221e+02     5.203210e-01
 * time: 7.555734872817993
    20     6.252215e+02     1.748976e-01
 * time: 7.8214709758758545
    21     6.252214e+02     3.233491e-02
 * time: 8.07071304321289
    22     6.252214e+02     3.132817e-02
 * time: 8.306206941604614
    23     6.252214e+02     7.487600e-02
 * time: 8.533815860748291
    24     6.252213e+02     8.901395e-02
 * time: 8.768064022064209
    25     6.252213e+02     7.637735e-02
 * time: 8.99588394165039
    26     6.252213e+02     5.953687e-02
 * time: 9.22209095954895
    27     6.252212e+02     5.053778e-02
 * time: 9.458745956420898
    28     6.252211e+02     9.222601e-02
 * time: 9.686851024627686
    29     6.252210e+02     1.014797e-01
 * time: 9.919749975204468
    30     6.252209e+02     6.355065e-02
 * time: 10.146324872970581
    31     6.252209e+02     1.650355e-02
 * time: 10.391231060028076
    32     6.252209e+02     1.160409e-02
 * time: 10.616590976715088
    33     6.252209e+02     9.804661e-03
 * time: 10.85079288482666
    34     6.252209e+02     1.127148e-02
 * time: 11.073854923248291
    35     6.252209e+02     9.136428e-03
 * time: 11.293308019638062
    36     6.252209e+02     3.817899e-03
 * time: 11.52663803100586
    37     6.252209e+02     3.746846e-03
 * time: 11.788095951080322

FittedPumasModel

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

Number of subjects:                             20

Observation records:         Active        Missing
    conc:                       160              0
    response:                   160              0
    Total:                      320              0

Number of parameters:      Constant      Optimized
                                  0             13

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                      NoXChange
Log-likelihood value:                    -625.2209

--------------------------
               Estimate
--------------------------
pop_Ka           0.94546
pop_CL           0.47437
pop_Vc          20.169
pop_R0         102.11
pop_turnover    10.219
pop_Imax         0.8705
pop_IC50         2.472
Ω_pk₁,₁          0.11826
Ω_pk₂,₂          0.08887
Ω_pk₃,₃          0.024739
Ω_pd₁,₁          0.028238
σ_prop_pk        0.098956
σ_add_pd         4.6475
--------------------------

coef(fit_result)

(pop_Ka = 0.9454609876410168,
 pop_CL = 0.474365974089112,
 pop_Vc = 20.16935944476136,
 pop_R0 = 102.10728122723762,
 pop_turnover = 10.219408933990612,
 pop_Imax = 0.8704951742242379,
 pop_IC50 = 2.4719877264988894,
 Ω_pk = [0.11825590901980063 0.0 0.0; 0.0 0.08887014391419974 0.0; 0.0 0.0 0.02473930839625612],
 Ω_pd = [0.02823775143489028;;],
 σ_prop_pk = 0.09895581112845298,
 σ_add_pd = 4.647498869249405,)

10 Model Diagnostics

insp = inspect(fit_result)

[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.

FittedPumasModelInspection

Likelihood approximation used for weighted residuals: FOCE

goodness_of_fit(insp; observations = [:response])

Figure 6: Goodness-of-fit for IDR Type I model

11 Summary

In this tutorial, we covered:

Turnover concept: Physiological variables maintained by production/removal balance
Four IDR types: Distinguished by drug action on k_{in} or k_{out}
Response direction: Types I/IV decrease response; Types II/III increase it
Turnover time: Governs speed of response and recovery
Pumas implementation: Complete models for all four types
Model selection: Based on drug mechanism of action

11.1 Key Takeaways

IDR models explain why effects outlast drug concentrations
The turnover half-life determines response dynamics
Choose model type based on mechanism (inhibition vs stimulation, production vs removal)
Effects are reversible - response returns to baseline after drug clears

11.2 Next Steps

If your drug shows:

Irreversible effects: Consider K-PD Models
Need to compare models: See Model Selection Tutorial

References

Chakraborty, A., W. Krzyzanski, and W. J. Jusko. 1999. “Mathematical Modeling of Circadian Cortisol Concentrations Using Indirect Response Models: Comparison of Several Methods.” Journal of Pharmacokinetics and Biopharmaceutics 27 (1): 23–43.

Dayneka, N. L., V. Garg, and W. J. Jusko. 1993. “Comparison of Four Basic Models of Indirect Pharmacodynamic Responses.” Journal of Pharmacokinetics and Biopharmaceutics 21 (4): 457–78.

Hamberg, A. K., M. L. Dahl, T. M. Barber, and M. Wadelius. 2007. “A PK-PD Model for Predicting the Impact of Age, CYP2C9, and VKORC1 Genotype on Individualization of Warfarin Therapy.” Clinical Pharmacology & Therapeutics 81 (4): 529–38.

Sharma, A., and W. J. Jusko. 1998. “Characteristics of Indirect Pharmacodynamic Models and Applications to Clinical Drug Responses.” British Journal of Clinical Pharmacology 45 (3): 229–39.

Reuse

CC BY-SA 4.0