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 ~ @. Normal(Cp, abs(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 ~ @. Normal(Cp, abs(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.060387e+02     7.932041e+01
 * time: 0.038537025451660156
     1     6.695555e+02     8.281340e+01
 * time: 1.7278509140014648
     2     6.486648e+02     3.484822e+01
 * time: 1.859673023223877
     3     6.407331e+02     3.158124e+01
 * time: 1.9816720485687256
     4     6.401870e+02     2.009610e+01
 * time: 2.0857229232788086
     5     6.390340e+02     2.041743e+01
 * time: 2.1951780319213867
     6     6.379018e+02     7.229470e+00
 * time: 2.300920009613037
     7     6.375381e+02     8.302646e+00
 * time: 2.404289960861206
     8     6.369474e+02     5.762145e+00
 * time: 2.6447970867156982
     9     6.364074e+02     6.339323e+00
 * time: 2.742219924926758
    10     6.361904e+02     6.503069e+00
 * time: 2.8403561115264893
    11     6.359350e+02     2.097937e+00
 * time: 2.9373910427093506
    12     6.358813e+02     1.098445e+00
 * time: 3.03645396232605
    13     6.358126e+02     1.083815e+00
 * time: 3.1406219005584717
    14     6.357031e+02     2.419836e+00
 * time: 3.2427101135253906
    15     6.356434e+02     1.338037e+00
 * time: 3.345870018005371
    16     6.356255e+02     4.260798e-01
 * time: 3.5046260356903076
    17     6.356197e+02     4.441574e-01
 * time: 3.5948519706726074
    18     6.356060e+02     1.051575e+00
 * time: 3.6858620643615723
    19     6.355902e+02     1.208733e+00
 * time: 3.777837038040161
    20     6.355771e+02     6.891400e-01
 * time: 3.87493896484375
    21     6.355715e+02     3.108437e-01
 * time: 3.977876901626587
    22     6.355677e+02     3.493174e-01
 * time: 4.077476978302002
    23     6.355610e+02     8.277929e-01
 * time: 4.199264049530029
    24     6.355499e+02     1.228623e+00
 * time: 4.296590089797974
    25     6.355365e+02     1.206323e+00
 * time: 4.391629934310913
    26     6.355278e+02     5.605031e-01
 * time: 4.489459037780762
    27     6.355254e+02     2.756764e-01
 * time: 4.5857579708099365
    28     6.355252e+02     3.877811e-02
 * time: 4.682050943374634
    29     6.355251e+02     3.410246e-02
 * time: 4.7877891063690186
    30     6.355249e+02     2.671407e-02
 * time: 4.877389907836914
    31     6.355249e+02     3.583602e-02
 * time: 4.967413902282715
    32     6.355248e+02     2.328782e-02
 * time: 5.055624961853027
    33     6.355248e+02     1.351482e-02
 * time: 5.1469550132751465
    34     6.355248e+02     1.432545e-02
 * time: 5.246841907501221
    35     6.355248e+02     2.232734e-02
 * time: 5.3352460861206055
    36     6.355247e+02     2.550356e-02
 * time: 5.42428994178772
    37     6.355247e+02     1.823888e-02
 * time: 5.522107124328613
    38     6.355247e+02     6.502648e-03
 * time: 5.610558032989502
    39     6.355247e+02     1.029964e-03
 * time: 5.699692010879517
    40     6.355247e+02     2.813378e-04
 * time: 5.79552698135376

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

-----------------------
            Estimate
-----------------------
pop_Ka        0.95244
pop_CL        0.52245
pop_Vc       19.685
pop_E0       53.82
pop_Emax    152.07
pop_EC50     10.849
Ω_pk₁,₁       0.033115
Ω_pk₂,₂       0.079942
Ω_pk₃,₃       0.073432
Ω_pd₁,₁       0.027538
Ω_pd₂,₂       0.016009
σ_prop_pk     0.095838
σ_add_pd      4.816
-----------------------

6.4 Inspect Results

# View estimated parameters
coef(fit_result)

(pop_Ka = 0.9524408424966573,
 pop_CL = 0.522446633138026,
 pop_Vc = 19.685141923230674,
 pop_E0 = 53.82020304789926,
 pop_Emax = 152.0729430302488,
 pop_EC50 = 10.848808510477031,
 Ω_pk = [0.033114819552790456 0.0 0.0; 0.0 0.07994216146356922 0.0; 0.0 0.0 0.0734321603045483],
 Ω_pd = [0.02753822847410206 0.0; 0.0 0.016008854415050303],
 σ_prop_pk = 0.09583769351496636,
 σ_add_pd = 4.816005081443679,)

# 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:                   -635.52473

---------------------------------------------------------------
            Estimate     SE          95.0% C.I.
---------------------------------------------------------------
pop_Ka        0.95244    0.064445    [  0.82613  ;   1.0788  ]
pop_CL        0.52245    0.033469    [  0.45685  ;   0.58804 ]
pop_Vc       19.685      1.2989      [ 17.139    ;  22.231   ]
pop_E0       53.82       4.768       [ 44.475    ;  63.165   ]
pop_Emax    152.07       61.724      [ 31.096    ; 273.05    ]
pop_EC50     10.849      7.8513      [ -4.5394   ;  26.237   ]
Ω_pk₁,₁       0.033115   0.025878    [ -0.017604 ;   0.083834]
Ω_pk₂,₂       0.079942   0.028441    [  0.024198 ;   0.13569 ]
Ω_pk₃,₃       0.073432   0.032399    [  0.0099322;   0.13693 ]
Ω_pd₁,₁       0.027538   0.01216     [  0.003706 ;   0.05137 ]
Ω_pd₂,₂       0.016009   0.022598    [ -0.028282 ;   0.0603  ]
σ_prop_pk     0.095838   0.0067094   [  0.082687 ;   0.10899 ]
σ_add_pd      4.816      0.25645     [  4.3134   ;   5.3186  ]
---------------------------------------------------------------

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.