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.

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

Figure 1: Direct response model: Effect is an immediate function of plasma concentration

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

Clinical 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

Potency 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

Figure 2: Effect of Hill coefficient (γ) on the concentration-effect curve

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 = u"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
end

PumasModel
  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, Effect

4.2 Understanding the Model Structure

Let’s break down the key components:

PK Component

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

PD Component

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}

Error Models

• 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 = 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 (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
end

PumasModel
  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, Response

5 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: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

Figure 3: Single-subject simulation showing parallel concentration and effect profiles

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

Figure 4: Population simulation showing between-subject variability in PK and PD

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),
)

Figure 5: Concentration-effect relationship showing no hysteresis (direct response)

No 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, effect

6.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.03715109825134277
     1     6.875973e+02     7.086690e+01
 * time: 1.7547359466552734
     2     6.730797e+02     2.770680e+01
 * time: 1.8857369422912598
     3     6.680278e+02     1.889767e+01
 * time: 2.007951021194458
     4     6.671751e+02     7.426803e+00
 * time: 2.7254719734191895
     5     6.663858e+02     1.148124e+01
 * time: 2.8232979774475098
     6     6.657989e+02     1.453872e+01
 * time: 2.9173269271850586
     7     6.644354e+02     5.981163e+00
 * time: 3.008197069168091
     8     6.639044e+02     1.549860e+00
 * time: 3.1047110557556152
     9     6.638476e+02     2.222896e+00
 * time: 3.244029998779297
    10     6.637499e+02     1.400233e+00
 * time: 3.3360860347747803
    11     6.635640e+02     1.374304e+00
 * time: 3.426598072052002
    12     6.634927e+02     8.691492e-01
 * time: 3.518440008163452
    13     6.634772e+02     4.672728e-01
 * time: 3.612821102142334
    14     6.634659e+02     4.958054e-01
 * time: 3.7422409057617188
    15     6.634489e+02     8.628975e-01
 * time: 3.8357880115509033
    16     6.634202e+02     1.077699e+00
 * time: 3.9275050163269043
    17     6.633958e+02     8.415645e-01
 * time: 4.018409967422485
    18     6.633845e+02     3.016369e-01
 * time: 4.123678922653198
    19     6.633811e+02     2.448260e-01
 * time: 4.216977119445801
    20     6.633781e+02     2.665923e-01
 * time: 4.309678077697754
    21     6.633729e+02     4.342333e-01
 * time: 4.4005279541015625
    22     6.633677e+02     5.116198e-01
 * time: 4.511608123779297
    23     6.633653e+02     2.433892e-01
 * time: 4.599771976470947
    24     6.633647e+02     1.213534e-01
 * time: 4.68770694732666
    25     6.633644e+02     7.003498e-02
 * time: 4.786942005157471
    26     6.633639e+02     1.218454e-01
 * time: 4.87444806098938
    27     6.633630e+02     2.684805e-01
 * time: 4.9620361328125
    28     6.633619e+02     2.732717e-01
 * time: 5.059546947479248
    29     6.633614e+02     1.187970e-01
 * time: 5.147257089614868
    30     6.633613e+02     3.184068e-02
 * time: 5.234940052032471
    31     6.633613e+02     2.863081e-02
 * time: 5.33238410949707
    32     6.633612e+02     4.296221e-02
 * time: 5.42051100730896
    33     6.633612e+02     4.155602e-02
 * time: 5.507122993469238
    34     6.633612e+02     1.877596e-02
 * time: 5.593228101730347
    35     6.633612e+02     3.084851e-03
 * time: 5.689391136169434
    36     6.633612e+02     3.079211e-03
 * time: 5.786827087402344
    37     6.633612e+02     3.077775e-03
 * time: 5.886718988418579
    38     6.633612e+02     3.077042e-03
 * time: 5.996537923812866
    39     6.633612e+02     3.076982e-03
 * time: 6.095556974411011
    40     6.633612e+02     3.076924e-03
 * time: 6.194478988647461
    41     6.633612e+02     3.076869e-03
 * time: 6.306222915649414
    42     6.633612e+02     3.076863e-03
 * time: 6.408260107040405
    43     6.633612e+02     3.076858e-03
 * time: 6.5106470584869385
    44     6.633612e+02     3.076857e-03
 * time: 6.620450973510742
    45     6.633612e+02     3.076857e-03
 * time: 6.718394994735718
    46     6.633612e+02     3.076856e-03
 * time: 6.81882905960083
    47     6.633612e+02     3.076856e-03
 * time: 6.931793928146362
    48     6.633612e+02     3.076856e-03
 * time: 7.038100004196167
    49     6.633612e+02     3.076856e-03
 * time: 7.144644021987915
    50     6.633612e+02     3.076856e-03
 * time: 7.273732900619507

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:              NoObjectiveChange
Log-likelihood value:                   -663.36118

-----------------------
            Estimate
-----------------------
pop_Ka        0.98228
pop_CL        0.49911
pop_Vc       20.097
pop_E0       52.319
pop_Emax    127.94
pop_EC50      8.2328
Ω_pk₁,₁       0.10792
Ω_pk₂,₂       0.079711
Ω_pk₃,₃       0.066919
Ω_pd₁,₁       0.034659
Ω_pd₂,₂       0.018256
σ_prop_pk     0.10323
σ_add_pd      5.266
-----------------------

6.4 Inspect Results

# View estimated parameters
coef(fit_result)

(pop_Ka = 0.982276167844883, pop_CL = 0.499112292572688, pop_Vc = 20.09706697620045, pop_E0 = 52.31877216550355, pop_Emax = 127.93520596859517, pop_EC50 = 8.23281529568625, Ω_pk = [0.1079207588701941 0.0 0.0; 0.0 0.07971053481870394 0.0; 0.0 0.0 0.06691876071568628], Ω_pd = [0.0346590654003878 0.0; 0.0 0.01825645373897743], σ_prop_pk = 0.10323021776261451, σ_add_pd = 5.266028562760617)

# 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:              NoObjectiveChange
Log-likelihood value:                   -663.36118

--------------------------------------------------------------
            Estimate     SE         95.0% C.I.
--------------------------------------------------------------
pop_Ka        0.98228    0.066785   [  0.85138  ;   1.1132  ]
pop_CL        0.49911    0.046677   [  0.40763  ;   0.5906  ]
pop_Vc       20.097      1.3206     [ 17.509    ;  22.685   ]
pop_E0       52.319      6.0667     [ 40.428    ;  64.209   ]
pop_Emax    127.94       31.672     [ 65.859    ; 190.01    ]
pop_EC50      8.2328     4.422      [ -0.43423  ;  16.9     ]
Ω_pk₁,₁       0.10792    0.051523   [  0.0069376;   0.2089  ]
Ω_pk₂,₂       0.079711   0.029916   [  0.021077 ;   0.13834 ]
Ω_pk₃,₃       0.066919   0.02633    [  0.015313 ;   0.11852 ]
Ω_pd₁,₁       0.034659   0.012566   [  0.010031 ;   0.059287]
Ω_pd₂,₂       0.018256   0.017476   [ -0.015997 ;   0.05251 ]
σ_prop_pk     0.10323    0.005469   [  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: FOCE

goodness_of_fit(insp; observations = [:conc])

Figure 6: Goodness-of-fit plot for concentration

goodness_of_fit(insp; observations = [:effect])

Figure 7: Goodness-of-fit plot for effect

7.2 Residual Diagnostics

wresiduals_vs_time(insp)

Figure 8: Weighted residuals over time

7.3 Visual Predictive Check

vpc_result = vpc(fit_result; observations = [:effect])

vpc_plot(vpc_result)

Figure 9: Visual Predictive Check for the effect observation

8 Summary

In this tutorial, we covered:

1. Conceptual foundation: Direct response models assume immediate equilibrium between concentration and effect
2. Mathematical framework: Emax, sigmoid Emax, and inhibitory formulations
3. Pumas implementation: Complete PK-PD models with BSV and RUV
4. Simulation: Single-subject and population simulations
5. Fitting: Parameter estimation using FOCE
6. 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, P. L. Gambus, C. Andresen, D. B. Goodale, and E. J. Youngs. 1999. “The Influence of Age on Propofol Pharmacodynamics.” Anesthesiology 90 (6): 1502–16.