flowchart LR
subgraph PK["Pharmacokinetics"]
Dose[Drug Dose] --> Absorption
Absorption --> Cp[Plasma<br/>Concentration Cp]
Cp --> Elimination
end
subgraph PD["Pharmacodynamics"]
Cp --> Effect["Effect<br/>E = E₀ + Emax·Cp/(EC₅₀+Cp)"]
end
Direct Response Pharmacodynamic Models
1 Learning Objectives
By the end of this tutorial, you will be able to:
- Understand when direct response models are appropriate
- Implement Emax and sigmoid Emax models in Pumas
- Simulate direct response PK-PD models
- Fit direct response models to data
- Perform model diagnostics and interpret results
- Distinguish between stimulatory and inhibitory formulations
2 Introduction
Direct response models represent the simplest form of pharmacodynamic modeling, where the drug effect is in immediate equilibrium with the plasma concentration. As the concentration rises, the effect rises; as the concentration falls, the effect falls - with no temporal delay between them.
2.1 When to Use Direct Response Models
Direct response models are appropriate when:
- The site of action is rapidly accessible from plasma
- Drug-receptor binding equilibrium is fast
- Signal transduction is nearly instantaneous
- No significant distribution delay exists
TipClinical Examples
Anesthetics: Propofol produces rapid onset and offset of sedation, closely tracking plasma concentrations (Schnider et al. 1999).
Bronchodilators: Short-acting beta-agonists like albuterol show rapid bronchodilation that parallels drug levels (Anderson 1993).
Antihypertensives: Some direct vasodilators produce immediate blood pressure changes proportional to concentration (Donnelly and Bhagat 1999).
2.2 When NOT to Use Direct Response Models
If you observe any of the following, consider alternative models:
- Hysteresis loops in effect vs. concentration plots → Effect compartment model
- Effect persists after drug concentration declines → Indirect response model
- Clockwise hysteresis → Tolerance or indirect mechanisms
3 Mathematical Framework
3.1 The Emax Model
The Emax model describes the concentration-effect relationship as a hyperbolic function (Holford and Sheiner 1981):
E = E_0 + \frac{E_{max} \cdot C}{EC_{50} + C}
| Parameter | Symbol | Description | Units |
|---|---|---|---|
| Baseline effect | E_0 | Effect without drug | Effect units |
| Maximum effect | E_{max} | Maximum drug-induced change | Effect units |
| Potency | EC_{50} | Concentration at 50% E_{max} | Concentration units |
| Concentration | C | Drug concentration | Concentration units |
NotePotency vs. Efficacy
EC_{50} measures potency: Lower values mean less drug is needed for effect. E_{max} measures efficacy: Higher values mean greater maximum effect. A drug can be highly potent (EC_{50} = 0.1 ng/mL) but have low efficacy (E_{max} = 10%).
3.2 The Sigmoid Emax (Hill) Model
When the concentration-effect relationship shows steepness different from hyperbolic, the Hill coefficient (\gamma) is added:
E = E_0 + \frac{E_{max} \cdot C^\gamma}{EC_{50}^\gamma + C^\gamma}
| Hill Coefficient (\gamma) | Interpretation |
|---|---|
| \gamma = 1 | Standard Emax (hyperbolic) |
| \gamma < 1 | Shallow curve, gradual response |
| \gamma > 1 | Steep curve, switch-like response |
| \gamma → \infty | All-or-none response |
3.3 Inhibitory Emax Model
For drugs that inhibit a response (e.g., enzyme inhibitors), the formulation is modified:
E = E_0 \cdot \left(1 - \frac{I_{max} \cdot C}{IC_{50} + C}\right)
Where:
- I_{max} = Maximum fractional inhibition (0 to 1)
- IC_{50} = Concentration producing 50% of maximum inhibition
4 Pumas Implementation
4.1 Stimulatory Direct Response Model
Let’s implement a complete PK-PD model with direct stimulatory response:
direct_response_stim = @model begin
@metadata begin
desc = "One-Compartment PK with Direct Emax PD Response"
timeu = "hr"
end
@param begin
# PK parameters
"""
Absorption rate constant (1/hr)
"""
pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
"""
Clearance (L/hr)
"""
pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
"""
Volume of distribution (L)
"""
pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)
# PD parameters
"""
Baseline response (no drug)
"""
pop_E0 ∈ RealDomain(; lower = 0, init = 50.0)
"""
Maximum drug-induced effect
"""
pop_Emax ∈ RealDomain(; lower = 0, init = 100.0)
"""
Concentration at 50% Emax (potency)
"""
pop_EC50 ∈ RealDomain(; lower = 0, init = 5.0)
# Between-subject variability
"""
Variance-covariance for PK parameters (CL, Vc, Ka)
"""
Ω_pk ∈ PDiagDomain(3)
"""
Variance-covariance for PD parameters (E0, Emax)
"""
Ω_pd ∈ PDiagDomain(2)
# Residual unexplained variability
"""
Proportional error for PK
"""
σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
"""
Additive error for PD
"""
σ_add_pd ∈ RealDomain(; lower = 0, init = 5.0)
end
@random begin
η_pk ~ MvNormal(Ω_pk)
η_pd ~ MvNormal(Ω_pd)
end
@pre begin
# Individual PK parameters
CL = pop_CL * exp(η_pk[1])
Vc = pop_Vc * exp(η_pk[2])
Ka = pop_Ka * exp(η_pk[3])
# Individual PD parameters
E0 = pop_E0 * exp(η_pd[1])
Emax = pop_Emax * exp(η_pd[2])
EC50 = pop_EC50 # Fixed at population level
end
@dynamics begin
Depot' = -Ka * Depot
Central' = Ka * Depot - (CL / Vc) * Central
end
@vars begin
# Plasma concentration (for plotting)
Cp = Central / Vc
# Direct Emax response (for plotting)
Effect = E0 + (Emax * Cp) / (EC50 + Cp)
end
@derived begin
# Observed values with residual error
conc ~ @. ProportionalNormal(Cp, σ_prop_pk)
effect ~ @. Normal(Effect, σ_add_pd)
end
endPumasModel
Parameters: pop_Ka, pop_CL, pop_Vc, pop_E0, pop_Emax, pop_EC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
Random effects: η_pk, η_pd
Covariates:
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived: conc, effect, Cp, Effect
Observed: conc, effect, Cp, Effect4.2 Understanding the Model Structure
Let’s break down the key components:
The pharmacokinetic model is a standard one-compartment model with first-order absorption:
\frac{d(Depot)}{dt} = -K_a \cdot Depot
\frac{d(Central)}{dt} = K_a \cdot Depot - \frac{CL}{V_c} \cdot Central
The pharmacodynamic component links concentration directly to effect:
C_p = \frac{Central}{V_c}
Effect = E_0 + \frac{E_{max} \cdot C_p}{EC_{50} + C_p}
- PK: Proportional error (larger concentrations → larger variability)
- PD: Additive error (constant variability across effect range)
4.3 Inhibitory Direct Response Model
For drugs that inhibit a response:
direct_response_inhib = @model begin
@metadata begin
desc = "One-Compartment PK with Direct Inhibitory Imax PD Response"
timeu = "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 (uninhibited)
"""
pop_R0 ∈ RealDomain(; lower = 0, init = 100.0)
"""
Maximum fractional inhibition (0 to 1)
"""
pop_Imax ∈ RealDomain(; lower = 0, upper = 1, init = 0.9)
"""
Concentration at 50% inhibition
"""
pop_IC50 ∈ RealDomain(; lower = 0, init = 5.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
CL = pop_CL * exp(η_pk[1])
Vc = pop_Vc * exp(η_pk[2])
Ka = pop_Ka * exp(η_pk[3])
R0 = pop_R0 * exp(η_pd[1])
Imax = pop_Imax
IC50 = pop_IC50
end
@dynamics begin
Depot' = -Ka * Depot
Central' = Ka * Depot - (CL / Vc) * Central
end
@vars begin
# Plasma concentration (for plotting)
Cp = Central / Vc
# Inhibitory Imax response (for plotting)
Inhibition = Imax * Cp / (IC50 + Cp)
Response = R0 * (1 - Inhibition)
end
@derived begin
conc ~ @. ProportionalNormal(Cp, σ_prop_pk)
response ~ @. Normal(Response, σ_add_pd)
end
endPumasModel
Parameters: pop_Ka, pop_CL, pop_Vc, pop_R0, pop_Imax, pop_IC50, Ω_pk, Ω_pd, σ_prop_pk, σ_add_pd
Random effects: η_pk, η_pd
Covariates:
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived: conc, response, Cp, Inhibition, Response
Observed: conc, response, Cp, Inhibition, Response5 Simulation Workflow
5.1 Single-Subject Simulation
Let’s simulate a single subject receiving an IV bolus dose:
# Define parameters
params_stim = (
pop_Ka = 1.0,
pop_CL = 0.5,
pop_Vc = 20.0,
pop_E0 = 50.0,
pop_Emax = 100.0,
pop_EC50 = 5.0,
Ω_pk = Diagonal([0.09, 0.09, 0.09]),
Ω_pd = Diagonal([0.04, 0.04]),
σ_prop_pk = 0.1,
σ_add_pd = 5.0,
)
# Define dosing regimen (100 mg oral dose)
dosing = DosageRegimen(100, time = 0, cmt = 1)
# Create subject
subject = Subject(id = 1, events = dosing)
# Simulate (without residual error for clean visualization)
sim_single = simobs(
direct_response_stim,
subject,
params_stim,
obstimes = 0:0.5:24,
simulate_error = false,
)SimulatedObservations
Simulated variables: conc, effect, Cp, Effect
Time: [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5 … 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0]5.2 Visualizing Concentration and Effect
# Convert simulation to DataFrame
sim_df = DataFrame(sim_single)
# Create PK plot (using vars: Cp)
pk_plt =
data(sim_df) *
mapping(:time => "Time (hr)", :Cp => "Concentration (mg/L)") *
visual(Lines, linewidth = 2, color = :steelblue)
# Create PD plot (using vars: Effect)
pd_plt =
data(sim_df) *
mapping(:time => "Time (hr)", :Effect => "Effect") *
visual(Lines, linewidth = 2, color = :coral)
fig = Figure(size = (900, 400))
draw!(fig[1, 1], pk_plt; axis = (title = "Plasma Concentration",))
draw!(fig[1, 2], pd_plt; axis = (title = "Pharmacodynamic Response",))
# Add baseline reference
hlines!(fig.content[2], [50], linestyle = :dash, color = :gray, label = "Baseline (E₀)")
fig
5.3 Population Simulation
To understand variability in response, let’s simulate a population:
# Create population of 20 subjects
population = map(1:20) do i
Subject(id = i, events = dosing)
end
# Simulate population
Random.seed!(123)
sim_pop = simobs(direct_response_stim, population, params_stim, obstimes = 0:0.5:24)Simulated population (Vector{<:Subject})
Simulated subjects: 20
Simulated variables: conc, effect, Cp, Effect# Convert population simulation to DataFrame
pop_df = DataFrame(sim_pop)
# Create PK plot with subjects (using vars: Cp)
pk_pop_plt =
data(pop_df) *
mapping(:time => "Time (hr)", :Cp => "Concentration (mg/L)", group = :id) *
visual(Lines, linewidth = 1, alpha = 0.5, color = :steelblue)
# Create PD plot with subjects (using vars: Effect)
pd_pop_plt =
data(pop_df) *
mapping(:time => "Time (hr)", :Effect => "Effect", group = :id) *
visual(Lines, linewidth = 1, alpha = 0.5, color = :coral)
fig = Figure(size = (900, 400))
draw!(fig[1, 1], pk_pop_plt; axis = (title = "Population PK Profiles",))
draw!(fig[1, 2], pd_pop_plt; axis = (title = "Population PD Profiles",))
fig
5.4 Concentration-Effect Relationship
A key diagnostic is plotting effect against concentration to verify the direct relationship:
# Use pop_df from previous code block
# Create concentration-effect plot
ce_plt =
data(pop_df) *
mapping(:conc => "Concentration (mg/L)", :effect => "Effect", group = :id) *
visual(Lines, linewidth = 1, alpha = 0.5, color = :steelblue)
# Add theoretical Emax curve
conc_range = 0:0.1:10
E0, Emax, EC50 = 50.0, 100.0, 5.0
theoretical_effect = E0 .+ (Emax .* conc_range) ./ (EC50 .+ conc_range)
theory_df = DataFrame(conc = collect(conc_range), effect = theoretical_effect)
theory_plt =
data(theory_df) *
mapping(:conc => "Concentration (mg/L)", :effect => "Effect") *
visual(Lines, linewidth = 3, color = :red, linestyle = :dash)
draw(
ce_plt + theory_plt;
axis = (title = "Concentration-Effect Relationship", width = 500, height = 400),
)
NoteNo Hysteresis in Direct Response
Notice that the concentration-effect curves overlay on the theoretical line without forming loops. This is the hallmark of a direct response model - effect follows concentration without delay. If you see loops (hysteresis), consider an effect compartment model.
6 Fitting Workflow
Let’s fit the direct response model to simulated data.
6.1 Generate Synthetic Data
# Simulate data with added noise
Random.seed!(456)
sim_data = simobs(
direct_response_stim,
population,
params_stim,
obstimes = [0.5, 1, 2, 4, 6, 8, 12, 24],
)
# Convert to DataFrame for fitting
df_sim = DataFrame(sim_data)
first(df_sim, 10)6.2 Create Population for Fitting
fit_pop = read_pumas(
df_sim;
id = :id,
time = :time,
amt = :amt,
evid = :evid,
cmt = :cmt,
observations = [:conc, :effect],
)Population
Subjects: 20
Observations: conc, effect6.3 Fit the Model
# Initial parameter estimates
init_params = (
pop_Ka = 0.8,
pop_CL = 0.4,
pop_Vc = 25.0,
pop_E0 = 45.0,
pop_Emax = 90.0,
pop_EC50 = 4.0,
Ω_pk = Diagonal([0.05, 0.05, 0.05]),
Ω_pd = Diagonal([0.02, 0.02]),
σ_prop_pk = 0.15,
σ_add_pd = 8.0,
)
# Fit using FOCE
fit_result = fit(direct_response_stim, 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.174606e+02 7.466666e+01 * time: 0.04540610313415527 1 6.878778e+02 6.530442e+01 * time: 9.21900200843811 2 6.746300e+02 3.072229e+01 * time: 9.278443098068237 3 6.690137e+02 2.302355e+01 * time: 9.331053972244263 4 6.674073e+02 8.907344e+00 * time: 9.380276918411255 5 6.668950e+02 1.200825e+01 * time: 9.420988082885742 6 6.660497e+02 1.180685e+01 * time: 10.554423093795776 7 6.653544e+02 8.831226e+00 * time: 10.588649988174438 8 6.645192e+02 3.904501e+00 * time: 10.622416019439697 9 6.638776e+02 2.887026e+00 * time: 10.656721115112305 10 6.638290e+02 3.813192e+00 * time: 10.690906047821045 11 6.637292e+02 1.282629e+00 * time: 10.724787950515747 12 6.636575e+02 9.869195e-01 * time: 10.758723020553589 13 6.635410e+02 2.348005e+00 * time: 10.794164896011353 14 6.634847e+02 1.448622e+00 * time: 10.832366943359375 15 6.634639e+02 4.212822e-01 * time: 10.870465993881226 16 6.634549e+02 4.109399e-01 * time: 10.909169912338257 17 6.634407e+02 9.601113e-01 * time: 10.946767091751099 18 6.634137e+02 1.366399e+00 * time: 10.984352111816406 19 6.633913e+02 1.004758e+00 * time: 11.021565914154053 20 6.633816e+02 3.618512e-01 * time: 11.05925703048706 21 6.633786e+02 2.473615e-01 * time: 11.096638917922974 22 6.633760e+02 3.516121e-01 * time: 11.13411808013916 23 6.633717e+02 6.064943e-01 * time: 11.171576023101807 24 6.633674e+02 5.655201e-01 * time: 11.291686058044434 25 6.633651e+02 2.358003e-01 * time: 11.324496030807495 26 6.633645e+02 6.336909e-02 * time: 11.361835956573486 27 6.633642e+02 9.164601e-02 * time: 11.394866943359375 28 6.633637e+02 1.936207e-01 * time: 11.427220106124878 29 6.633629e+02 2.664015e-01 * time: 11.45953106880188 30 6.633619e+02 2.341386e-01 * time: 11.492435932159424 31 6.633614e+02 1.041297e-01 * time: 11.526124000549316 32 6.633612e+02 2.508651e-02 * time: 11.562709093093872 33 6.633612e+02 2.197022e-02 * time: 11.599008083343506 34 6.633612e+02 1.797574e-02 * time: 11.635973930358887 35 6.633612e+02 1.547882e-02 * time: 11.671641111373901 36 6.633612e+02 1.160593e-02 * time: 11.707424879074097 37 6.633612e+02 2.969454e-03 * time: 11.778614044189453 38 6.633612e+02 1.082502e-03 * time: 11.809149980545044 39 6.633612e+02 1.585648e-04 * time: 11.838582992553711
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 20
Observation records: Active Missing
conc: 160 0
effect: 160 0
Total: 320 0
Number of parameters: Constant Optimized
0 13
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -663.36119
-----------------------
Estimate
-----------------------
pop_Ka 0.98228
pop_CL 0.49911
pop_Vc 20.097
pop_E0 52.317
pop_Emax 127.92
pop_EC50 8.2315
Ω_pk₁,₁ 0.10792
Ω_pk₂,₂ 0.079712
Ω_pk₃,₃ 0.066925
Ω_pd₁,₁ 0.034658
Ω_pd₂,₂ 0.018256
σ_prop_pk 0.10323
σ_add_pd 5.266
-----------------------6.4 Inspect Results
# View estimated parameters
coef(fit_result)(pop_Ka = 0.9822770152375645, pop_CL = 0.4991110615321096, pop_Vc = 20.096961450169967, pop_E0 = 52.317447853505186, pop_Emax = 127.92449497890178, pop_EC50 = 8.231470287370728, Ω_pk = [0.10792019145465229 0.0 0.0; 0.0 0.07971248904729981 0.0; 0.0 0.0 0.06692511900209187], Ω_pd = [0.03465771778074073 0.0; 0.0 0.018255542319875622], σ_prop_pk = 0.10322966487465185, σ_add_pd = 5.266031963843887)# View confidence intervals
infer(fit_result)[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: Matrix exponential
Number of subjects: 20
Observation records: Active Missing
conc: 160 0
effect: 160 0
Total: 320 0
Number of parameters: Constant Optimized
0 13
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -663.36119
--------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------
pop_Ka 0.98228 0.066786 [ 0.85138 ; 1.1132 ]
pop_CL 0.49911 0.046677 [ 0.40763 ; 0.5906 ]
pop_Vc 20.097 1.3205 [ 17.509 ; 22.685 ]
pop_E0 52.317 6.066 [ 40.428 ; 64.207 ]
pop_Emax 127.92 31.655 [ 65.881 ; 189.97 ]
pop_EC50 8.2315 4.4198 [ -0.43125 ; 16.894 ]
Ω_pk₁,₁ 0.10792 0.051522 [ 0.006938; 0.2089 ]
Ω_pk₂,₂ 0.079712 0.029917 [ 0.021077; 0.13835 ]
Ω_pk₃,₃ 0.066925 0.026334 [ 0.015312; 0.11854 ]
Ω_pd₁,₁ 0.034658 0.012565 [ 0.010032; 0.059284]
Ω_pd₂,₂ 0.018256 0.017474 [ -0.015994; 0.052505]
σ_prop_pk 0.10323 0.0054689 [ 0.092511; 0.11395 ]
σ_add_pd 5.266 0.37852 [ 4.5241 ; 6.0079 ]
--------------------------------------------------------------7 Model Diagnostics
7.1 Goodness-of-Fit Plots
# Inspect the fit
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: FOCEgoodness_of_fit(insp; observations = [:conc])
goodness_of_fit(insp; observations = [:effect])
7.2 Residual Diagnostics
wresiduals_vs_time(insp)
7.3 Visual Predictive Check
vpc_result = vpc(fit_result; observations = [:effect])vpc_plot(vpc_result)
8 Summary
In this tutorial, we covered:
- Conceptual foundation: Direct response models assume immediate equilibrium between concentration and effect
- Mathematical framework: Emax, sigmoid Emax, and inhibitory formulations
- Pumas implementation: Complete PK-PD models with BSV and RUV
- Simulation: Single-subject and population simulations
- Fitting: Parameter estimation using FOCE
- Diagnostics: GOF plots, residuals, and VPC
8.1 Key Takeaways
- Direct response is the simplest PD model but only appropriate when effect equilibrates rapidly
- The EC_{50} parameter characterizes potency, E_{max} characterizes efficacy
- Absence of hysteresis is the diagnostic signature of direct response
- Both stimulatory and inhibitory formulations follow the same principles
8.2 Next Steps
If your data shows:
- Hysteresis: Proceed to Effect Compartment Models
- Effect persists after drug clears: Proceed to Indirect Response Models
- Model selection challenges: See Model Selection Tutorial
References
Anderson, G. P. 1993. “Formoterol: Pharmacology, Molecular Basis of Agonism, and Mechanism of Long Duration of a Highly Potent and Selective Beta2-Adrenoceptor Agonist Bronchodilator.” Life Sciences 52 (26): 2145–60.
Donnelly, R., and K. Bhagat. 1999. “Direct Vascular Effects of Calcium Antagonists.” Clinical and Experimental Hypertension 21 (5-6): 969–85.
Holford, N. H. G., and L. B. Sheiner. 1981. “Understanding the Dose-Effect Relationship: Clinical Application of Pharmacokinetic-Pharmacodynamic Models.” Clinical Pharmacokinetics 6 (6): 429–53.
Schnider, T. W., C. F. Minto, S. L. Shafer, et al. 1999. “The Influence of Age on Propofol Pharmacodynamics.” Anesthesiology 90 (6): 1502–16.