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 = 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
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 = 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
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:24.05.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.04110598564147949 1 6.875973e+02 7.086690e+01 * time: 2.658334970474243 2 6.730797e+02 2.770680e+01 * time: 3.45778489112854 3 6.680278e+02 1.889767e+01 * time: 3.5883820056915283 4 6.671751e+02 7.426803e+00 * time: 3.70209002494812 5 6.663858e+02 1.148124e+01 * time: 3.8069798946380615 6 6.657989e+02 1.453872e+01 * time: 3.9130640029907227 7 6.644354e+02 5.981163e+00 * time: 4.076335906982422 8 6.639044e+02 1.549860e+00 * time: 4.1758129596710205 9 6.638476e+02 2.222896e+00 * time: 4.27664589881897 10 6.637499e+02 1.400233e+00 * time: 4.37617301940918 11 6.635640e+02 1.374304e+00 * time: 4.47855806350708 12 6.634927e+02 8.691492e-01 * time: 4.6003758907318115 13 6.634772e+02 4.672728e-01 * time: 4.697490930557251 14 6.634659e+02 4.958054e-01 * time: 4.7942070960998535 15 6.634489e+02 8.628975e-01 * time: 4.887834072113037 16 6.634202e+02 1.077699e+00 * time: 4.99358606338501 17 6.633958e+02 8.415645e-01 * time: 5.086967945098877 18 6.633845e+02 3.016369e-01 * time: 5.181713104248047 19 6.633811e+02 2.448260e-01 * time: 5.288470029830933 20 6.633781e+02 2.665923e-01 * time: 5.384721994400024 21 6.633729e+02 4.342333e-01 * time: 5.479794025421143 22 6.633677e+02 5.116198e-01 * time: 5.574812889099121 23 6.633653e+02 2.433892e-01 * time: 5.678276062011719 24 6.633647e+02 1.213534e-01 * time: 5.77349591255188 25 6.633644e+02 7.003498e-02 * time: 5.865998983383179 26 6.633639e+02 1.218454e-01 * time: 5.969944000244141 27 6.633630e+02 2.684805e-01 * time: 6.064276933670044 28 6.633619e+02 2.732717e-01 * time: 6.156876087188721 29 6.633614e+02 1.187970e-01 * time: 6.26275110244751 30 6.633613e+02 3.184068e-02 * time: 6.355009078979492 31 6.633613e+02 2.863081e-02 * time: 6.449107885360718 32 6.633612e+02 4.296221e-02 * time: 6.55264687538147 33 6.633612e+02 4.155602e-02 * time: 6.650842905044556 34 6.633612e+02 1.877596e-02 * time: 6.747485876083374 35 6.633612e+02 3.084851e-03 * time: 6.845194101333618 36 6.633612e+02 3.079211e-03 * time: 6.963658094406128 37 6.633612e+02 3.077775e-03 * time: 7.076896905899048 38 6.633612e+02 3.077042e-03 * time: 7.189709901809692 39 6.633612e+02 3.076982e-03 * time: 7.3159520626068115 40 6.633612e+02 3.076925e-03 * time: 7.430368900299072 41 6.633612e+02 3.076869e-03 * time: 7.542666912078857 42 6.633612e+02 3.076863e-03 * time: 7.669591903686523 43 6.633612e+02 3.076858e-03 * time: 7.787549018859863 44 6.633612e+02 3.076857e-03 * time: 7.918094873428345 45 6.633612e+02 3.076857e-03 * time: 8.03563904762268 46 6.633612e+02 3.076856e-03 * time: 8.1545569896698 47 6.633612e+02 3.076856e-03 * time: 8.288832902908325 48 6.633612e+02 3.076856e-03 * time: 8.409761905670166 49 6.633612e+02 3.076856e-03 * time: 8.530729055404663 50 6.633612e+02 3.076856e-03 * time: 8.663424968719482 51 6.633612e+02 3.076856e-03 * time: 8.791853904724121 52 6.633612e+02 3.076856e-03 * time: 8.945453882217407 53 6.633612e+02 3.076856e-03 * time: 9.088401079177856
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.9822761678449141, pop_CL = 0.4991122925727062, pop_Vc = 20.097066976201404, pop_E0 = 52.31877216550632, pop_Emax = 127.93520596865199, pop_EC50 = 8.23281529569161, Ω_pk = [0.10792075886985882 0.0 0.0; 0.0 0.07971053481872022 0.0; 0.0 0.0 0.0669187607156609], Ω_pd = [0.03465906540038306 0.0; 0.0 0.018256453738992706], σ_prop_pk = 0.10323021776261089, σ_add_pd = 5.266028562760324)# 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: 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, 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.