MI-210 Essentials of Population PKPD M&S — Pumas Edition
This chapter covers the foundational concepts for population PK-PD modeling:
Population pharmacokinetic-pharmacodynamic analysis aims to:
Pharmacometric analysis rests on two stacked pillars. The lower pillar is technical expertise — statistics, mathematics, and programming — accounting for roughly 20–30% of the analytical effort. The upper pillar is strategic expertise — scientific judgment, question formulation, and decision-making — which accounts for 70–80%. Without the technical foundation the strategic layer has no basis, but technical skill alone is insufficient.
The principle of parsimony is central: the goal is the simplest model that adequately describes the central tendency and variability of the data. Not every analysis requires a full NLME model; some questions are best answered by NCA. The complexity of the model should be commensurate with the complexity of the question.
As John Tukey observed: “It is better to find an approximate answer to the right question than to find the right answer to an approximate question.” Investing time in formulating the right question — grounded in knowledge of the drug, disease, and trial — is more valuable than technical precision applied to a poorly defined problem.
Goal: To understand factors leading to variability in pharmacokinetic and pharmacodynamic response for appropriate drug use.
“Drug Development = Model-Building”
A multiple cohort ascending single dose study of MI2005A. One endpoint of interest is the effect on cardiac repolarization (delta QTc vs. plasma concentration).
dqtc_df = CSV.read(joinpath(data_dir, "dqtc.csv"), DataFrame)
rename!(dqtc_df, "SID" => :id, "CONC" => :conc, "dqtc" => :dqtc, "DOSE" => :dose, "TIME" => :time)
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1],
xlabel = "Concentration (ng/mL)",
ylabel = "Delta QTc (msec)",
)
scatter!(ax, dqtc_df.conc, dqtc_df.dqtc, color = :steelblue, markersize = 6, alpha = 0.6)
fig
Fitting a Model to Data:
A parameter (e.g., \(k_e\), \(V\)) is a numerical quantity estimated from observed data that defines the relationship between independent and dependent variables within a given model. Once estimated, it is treated as fixed within that model context. A constant is a value held fixed from prior knowledge — for example, hepatic blood flow (\(\approx 1.5\) L/hr/kg) is treated as a physiological constant in many PK models, even though it was itself originally estimated as a parameter from experimental data. This distinction is frequently conflated in the literature.
All data contain error, and this extends beyond model residuals:
Identifying and characterizing these error sources before fitting any model is a critical and often underappreciated step in the analysis workflow.
Let observed data \(Y\) be described with parameter \(\theta\). The probability \(P(Y)\) can be modeled as a function of \(\theta\).
The maximum likelihood estimate (MLE) \(\hat{\theta}\) is the value of \(\theta\) that maximizes \(L(Y|\theta)\).
For normally-distributed measurements:
\[y_j \sim N(f(x_j, \theta), \sigma^2)\]
The ML objective function (equivalent to \(-2\log L\) minus constants):
\[OF_{ML}(x_j, \theta, \sigma^2) = \sum_{j=1}^{n} \left[ \frac{(y_j - f(x_j, \theta))^2}{\sigma^2} + \log \sigma^2 \right]\]
Ordinary Least Squares (OLS):
\[OF_{OLS} = \sum_{j=1}^{n} (obs_j - pred_j)^2\]
Weighted Least Squares (WLS) (where \(W\) is typically \(1/obs\)):
\[OF_{WLS} = \sum_{j=1}^{n} \left[ (obs_j - pred_j)^2 \cdot W_j \right]\]
Extended Least Squares (ELS) (Maximum Likelihood):
\[OF_{ELS} = \sum_{j=1}^{n} \left[ \frac{(obs_j - pred_j)^2}{\sigma_j^2} + \log \sigma_j^2 \right]\]
The book uses an Excel workbook to demonstrate nonlinear regression on a 1-compartment IV bolus PK model. Here we reproduce this in Julia/Pumas.
Model: \(C_t = \frac{Dose}{V} \cdot e^{-\frac{CL}{V} \cdot t}\)
True values: CL = 3.0 L/hr, V = 20.0 L
The one-compartment model treats the entire body as a single, well-mixed space of volume \(V\). Following an IV bolus, the drug is assumed to distribute instantaneously throughout this space, yielding an initial concentration \(C_0 = \text{Dose}/V\). Thereafter, concentration declines monoexponentially according to first-order elimination:
\[C(t) = \frac{\text{Dose}}{V} \cdot e^{-k_e \cdot t}\]
where the elimination rate constant \(k_e = CL/V\) is obtained from the slope of the log-linear concentration–time profile. Clearance \(CL\) governs the rate of irreversible drug removal (units: volume/time), while \(V\) relates the amount of drug in the body to the measured plasma concentration. On a semi-logarithmic plot, a straight line confirms monoexponential decline and supports the one-compartment assumption. Deviations from linearity (e.g., a biphasic decline) indicate distribution into peripheral tissues and motivate multi-compartment models.
# Simulate data from 1-cpt IV bolus (as in the Excel workbook)
# True parameters: CL=3, V=20, Dose=1000, additive sigma=0.9
true_CL = 3.0
true_V = 20.0
dose = 1000.0
times = [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 18, 19]
true_conc = @. dose / true_V * exp(-true_CL / true_V * times)
# Add noise (from book's dataset — heteroscedastic, proportional to prediction)
using Random
Random.seed!(210)
obs_conc = true_conc .+ randn(length(times)) .* (0.9 .* sqrt.(true_conc))
pk_data = DataFrame(
time = Float64.(times),
obs = obs_conc,
pred_true = true_conc,
)| Row | time | obs | pred_true |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 1.0 | 47.1589 | 43.0354 |
| 2 | 2.0 | 34.7978 | 37.0409 |
| 3 | 3.0 | 24.1227 | 31.8814 |
| 4 | 4.0 | 26.7429 | 27.4406 |
| 5 | 5.0 | 17.734 | 23.6183 |
| 6 | 6.0 | 17.3789 | 20.3285 |
| 7 | 7.0 | 16.0105 | 17.4969 |
| 8 | 8.0 | 9.12062 | 15.0597 |
| 9 | 10.0 | 11.3269 | 11.1565 |
| 10 | 12.0 | 11.0079 | 8.26494 |
| 11 | 14.0 | 8.23576 | 6.12282 |
| 12 | 16.0 | 2.58291 | 4.5359 |
| 13 | 18.0 | 0.99968 | 3.36028 |
| 14 | 19.0 | -1.29063 | 2.89222 |
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Time (hr)", ylabel = "Concentration (mg/L)")
scatter!(ax, pk_data.time, pk_data.obs, label = "Observed", color = :blue)
lines!(ax, pk_data.time, pk_data.pred_true, label = "True", color = :magenta, linewidth = 2)
axislegend(ax, position = :rt)
fig
We fit the 1-cpt IV bolus model using Pumas to find the MLE of CL and V.
# Create a Pumas-compatible dataset
# Single subject, IV bolus at time 0
ols_df = DataFrame(
id = fill(1, length(times) + 1),
time = vcat(0.0, Float64.(times)),
dv = vcat(missing, obs_conc),
amt = vcat(dose, fill(0.0, length(times))),
evid = vcat(1, fill(0, length(times))),
cmt = fill(1, length(times) + 1),
)
ols_pop = read_pumas(ols_df)Population
Subjects: 1
Observations: dv# 1-compartment IV bolus model with additive error
pk_model_add = @model begin
@param begin
tvCL ∈ RealDomain(lower = 0.0, init = 5.0)
tvV ∈ RealDomain(lower = 0.0, init = 30.0)
σ ∈ RealDomain(lower = 0.0, init = 1.0)
end
@pre begin
CL = tvCL
Vc = tvV
end
@dynamics Central1
@derived begin
cp := @. Central / Vc
dv ~ @. Normal(cp, σ)
end
endPumasModel
Parameters: tvCL, tvV, σ
Random effects:
Covariates:
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: dv
Observed: dvols_params = (tvCL = 5.0, tvV = 30.0, σ = 1.0)
ols_fit = fit(pk_model_add, ols_pop, ols_params, Pumas.NaivePooled())
ols_fit[ 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 3.902441e+02 7.407580e+02 * time: 0.013134002685546875 1 1.838690e+02 4.303081e+02 * time: 0.9699008464813232 2 4.315205e+01 3.654593e+01 * time: 0.9701130390167236 3 4.130316e+01 2.856668e+01 * time: 0.9702048301696777 4 3.754513e+01 6.837781e+00 * time: 0.97027587890625 5 3.721918e+01 6.838372e+00 * time: 0.9703469276428223 6 3.707849e+01 6.475951e+00 * time: 0.9704129695892334 7 3.682665e+01 5.802590e+00 * time: 0.970477819442749 8 3.646344e+01 7.886149e+00 * time: 0.9705429077148438 9 3.611610e+01 8.065761e+00 * time: 0.9706118106842041 10 3.600781e+01 4.828717e+00 * time: 0.9706759452819824 11 3.594982e+01 4.184895e+00 * time: 0.9707398414611816 12 3.588530e+01 1.296320e+00 * time: 0.9708058834075928 13 3.587922e+01 3.708834e-01 * time: 0.9708728790283203 14 3.587873e+01 1.223761e-02 * time: 0.9709658622741699 15 3.587872e+01 9.893961e-04 * time: 0.9710628986358643
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 1
Observation records: Active Missing
dv: 14 0
Total: 14 0
Number of parameters: Constant Optimized
0 3
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -35.878724
----------------
Estimate
----------------
tvCL 3.5122
tvV 19.451
σ 3.1388
----------------pk_model_prop = @model begin
@param begin
tvCL ∈ RealDomain(lower = 0.0, init = 5.0)
tvV ∈ RealDomain(lower = 0.0, init = 30.0)
σprop ∈ RealDomain(lower = 0.0, init = 0.2)
end
@pre begin
CL = tvCL
Vc = tvV
end
@dynamics Central1
@derived begin
cp := @. Central / Vc
dv ~ @. Normal(cp, abs(cp) * σprop)
end
end
prop_params = (tvCL = 5.0, tvV = 30.0, σprop = 0.2)
prop_fit = fit(pk_model_prop, ols_pop, prop_params, Pumas.NaivePooled())
prop_fit[ 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 1.497760e+02 6.489612e+02 * time: 1.8835067749023438e-5 1 6.583558e+01 5.023471e+01 * time: 0.26428890228271484 2 6.042026e+01 3.609868e+01 * time: 0.26444292068481445 3 5.405475e+01 1.658516e+01 * time: 0.26451992988586426 4 5.202441e+01 9.374855e+00 * time: 0.26459383964538574 5 5.106977e+01 9.911111e+00 * time: 0.264664888381958 6 5.059487e+01 1.036664e+01 * time: 0.2647368907928467 7 5.006839e+01 1.102774e+01 * time: 0.26480698585510254 8 4.895126e+01 1.245972e+01 * time: 0.2648749351501465 9 4.844996e+01 1.287806e+01 * time: 0.2649509906768799 10 4.721462e+01 1.188565e+01 * time: 0.2650458812713623 11 4.681847e+01 8.954011e+00 * time: 0.26515698432922363 12 4.652645e+01 6.250103e+00 * time: 0.26524996757507324 13 4.591016e+01 6.234626e+00 * time: 0.26540684700012207 14 4.483006e+01 3.945992e+00 * time: 0.2654998302459717 15 4.479167e+01 5.529813e+00 * time: 0.2655768394470215 16 4.455926e+01 2.417749e+00 * time: 0.26567602157592773 17 4.453781e+01 2.781958e+00 * time: 0.26576995849609375 18 4.450546e+01 3.241736e+00 * time: 0.2658510208129883 19 4.445377e+01 3.518406e+00 * time: 0.26592087745666504 20 4.432765e+01 3.360853e+00 * time: 0.26599693298339844 21 4.420763e+01 1.998236e+00 * time: 0.2660658359527588 22 4.413724e+01 7.509550e-01 * time: 0.26613593101501465 23 4.412525e+01 1.131628e-01 * time: 0.2662038803100586 24 4.412495e+01 3.887471e-03 * time: 0.26627182960510254 25 4.412495e+01 1.983269e-04 * time: 0.2663400173187256
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 1
Observation records: Active Missing
dv: 14 0
Total: 14 0
Number of parameters: Constant Optimized
0 3
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -44.124954
-----------------
Estimate
-----------------
tvCL 3.8249
tvV 28.052
σprop 0.5361
-----------------Standard diagnostic plots for evaluating model fit:
insp = inspect(ols_fit)
insp_df = DataFrame(insp)
obs_rows = @rsubset(insp_df, !ismissing(:dv))
fig = Figure(size = (800, 600))
# PRED vs OBS
ax1 = Axis(fig[1, 1], xlabel = "Observed", ylabel = "Predicted", title = "PRED vs OBS")
scatter!(ax1, obs_rows.dv, obs_rows.dv_pred, color = :navy, markersize = 6)
ablines!(ax1, 0, 1, color = :black)
# RES vs PRED
ax2 = Axis(fig[1, 2], xlabel = "Predicted", ylabel = "Residuals", title = "RES vs PRED")
res = obs_rows.dv .- obs_rows.dv_pred
scatter!(ax2, obs_rows.dv_pred, res, color = :navy, markersize = 6)
hlines!(ax2, 0, color = :black)
# RES vs TIME
ax3 = Axis(fig[2, 1], xlabel = "Time", ylabel = "Residuals", title = "RES vs TIME")
scatter!(ax3, obs_rows.time, res, color = :navy, markersize = 6)
hlines!(ax3, 0, color = :black)
# WRES vs PRED
ax4 = Axis(fig[2, 2], xlabel = "Predicted", ylabel = "Weighted Residuals", title = "WRES vs PRED")
wres = res ./ coef(ols_fit).σ
scatter!(ax4, obs_rows.dv_pred, wres, color = :navy, markersize = 6)
hlines!(ax4, 0, color = :black)
fig[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
The book compares parameter estimates across estimation methods:
| Parameter | LR | OLS | WLS | ELS | True |
|---|---|---|---|---|---|
| CL (L/hr) | 3.18 | 3.11 | 3.28 | 3.09 | 3.00 |
| V (L) | 25.52 | 24.82 | 25.14 | 25.30 | 20.00 |
When data are heteroscedastic (variance proportional to the mean), ELS/WLS give better estimates than OLS. In Pumas, this is handled by choosing the appropriate error model (proportional vs. additive) in the @derived block.
Using the Phase 1 MI2005A delta QTc data, fit a linear model (naive pooled approach) and compare OLS vs. ELS.
# Keep original subject IDs — NaivePooled estimation ignores grouping structure
# Use row index within each subject as the time axis (CONC is the covariate)
dqtc_pool = DataFrame(
id = dqtc_df.id,
conc = Float64.(dqtc_df.conc),
dv = Float64.(dqtc_df.dqtc),
)
dqtc_pool.time = Vector{Float64}(undef, nrow(dqtc_pool))
for gdf in groupby(dqtc_pool, :id)
gdf.time .= Float64.(1:nrow(gdf))
end
dqtc_pop = read_pumas(
dqtc_pool;
observations = [:dv],
covariates = [:conc],
event_data = false,
)Population
Subjects: 64
Covariates: conc
Observations: dv# Linear model: dQTc = intercept + slope * CONC
linear_dqtc = @model begin
@param begin
tvINT ∈ RealDomain(init = 1.0)
tvSLP ∈ RealDomain(init = 0.02)
σ_add ∈ RealDomain(lower = 0.0, init = 5.0)
end
@covariates conc
@pre begin
INT = tvINT
SLP = tvSLP
end
@derived begin
mu := @. INT + SLP * conc
dv ~ @. Normal(mu, σ_add)
end
end
lin_params = (tvINT = 1.0, tvSLP = 0.02, σ_add = 5.0)
lin_fit = fit(linear_dqtc, dqtc_pop, lin_params, Pumas.NaivePooled())
lin_fit[ 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 2.177739e+03 1.418514e+04 * time: 2.002716064453125e-5 1 2.145965e+03 6.258659e+02 * time: 0.405163049697876 2 2.118681e+03 7.104756e+02 * time: 0.4053988456726074 3 1.915684e+03 2.055385e+03 * time: 0.4056978225708008 4 1.856796e+03 2.183562e+03 * time: 0.4059309959411621 5 1.856667e+03 8.088581e+03 * time: 0.40617799758911133 6 1.853771e+03 2.659683e+03 * time: 0.40639686584472656 7 1.851798e+03 1.320097e+03 * time: 0.40663790702819824 8 1.847591e+03 1.163686e+02 * time: 0.40679287910461426 9 1.814822e+03 6.855837e+03 * time: 0.406965970993042 10 1.791657e+03 8.789379e+03 * time: 0.4071519374847412 11 1.785043e+03 3.339229e+03 * time: 0.4072990417480469 12 1.783504e+03 1.702782e+02 * time: 0.4074678421020508 13 1.783244e+03 1.914953e+02 * time: 0.4076099395751953 14 1.783242e+03 2.541537e+00 * time: 0.40775585174560547 15 1.783242e+03 1.408389e-01 * time: 0.4079718589782715 16 1.783242e+03 3.340979e-03 * time: 0.4081108570098877 17 1.783242e+03 6.584688e-05 * time: 0.40825581550598145
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 64
Observation records: Active Missing
dv: 812 0
Total: 812 0
Number of parameters: Constant Optimized
0 3
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1783.2424
------------------
Estimate
------------------
tvINT -0.088942
tvSLP 0.016926
σ_add 2.1753
------------------lin_coefs = coef(lin_fit)
conc_range = range(0, maximum(dqtc_df.conc), length = 100)
pred_dqtc = lin_coefs.tvINT .+ lin_coefs.tvSLP .* conc_range
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Concentration (ng/mL)", ylabel = "Delta QTc (msec)")
scatter!(ax, dqtc_df.conc, dqtc_df.dqtc, color = :steelblue, markersize = 5, alpha = 0.5)
lines!(ax, collect(conc_range), pred_dqtc, color = :red, linewidth = 2, label = "Linear Fit")
axislegend(ax, position = :lt)
fig
Fit a linear and Emax model to the pooled Phase 1 MI2005A AST data.
ast_df = CSV.read(joinpath(data_dir, "ast.csv"), DataFrame)
rename!(ast_df, "SID" => :id, "CONC" => :conc, "AST" => :ast, "DOSE" => :dose, "TIME" => :time)
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Concentration (ng/mL)", ylabel = "AST")
scatter!(ax, ast_df.conc, ast_df.ast, color = :steelblue, markersize = 5, alpha = 0.5)
fig
# Keep original subject IDs, use row index as time
ast_pool = DataFrame(
id = ast_df.id,
conc = Float64.(ast_df.conc),
dv = Float64.(ast_df.ast),
)
ast_pool.time = Vector{Float64}(undef, nrow(ast_pool))
for gdf in groupby(ast_pool, :id)
gdf.time .= Float64.(1:nrow(gdf))
end
ast_pop = read_pumas(
ast_pool;
observations = [:dv],
covariates = [:conc],
event_data = false,
)
# Emax model: AST = E0 + Emax*CONC/(EC50 + CONC)
emax_ast = @model begin
@param begin
tvE0 ∈ RealDomain(init = 15.0)
tvEmax ∈ RealDomain(init = 50.0)
tvEC50 ∈ RealDomain(lower = 0.0, init = 200.0)
σ_add ∈ RealDomain(lower = 0.0, init = 5.0)
end
@covariates conc
@pre begin
E0 = tvE0
Emax = tvEmax
EC50 = tvEC50
end
@derived begin
mu := @. E0 + Emax * conc / (EC50 + conc)
dv ~ @. Normal(mu, σ_add)
end
end
emax_params = (tvE0 = 15.0, tvEmax = 50.0, tvEC50 = 200.0, σ_add = 5.0)
emax_fit = fit(emax_ast, ast_pop, emax_params, Pumas.NaivePooled())
emax_fit[ 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 9.558379e+02 4.883388e+02 * time: 2.4080276489257812e-5 1 8.734043e+02 1.465283e+02 * time: 0.5386660099029541 2 8.449657e+02 9.520802e+01 * time: 0.5389299392700195 3 8.380284e+02 9.614958e+01 * time: 0.5391139984130859 4 8.306316e+02 1.819355e+01 * time: 0.5392608642578125 5 8.298644e+02 1.097222e+01 * time: 0.5393779277801514 6 8.294326e+02 7.037617e+00 * time: 0.5394999980926514 7 8.292394e+02 6.448246e+00 * time: 0.5396370887756348 8 8.291609e+02 1.623736e+00 * time: 0.5397610664367676 9 8.291570e+02 4.450326e-01 * time: 0.5398809909820557 10 8.291558e+02 4.445201e-01 * time: 0.5399999618530273 11 8.291504e+02 1.461467e+00 * time: 0.5401198863983154 12 8.291394e+02 2.878939e+00 * time: 0.5402369499206543 13 8.291117e+02 4.926721e+00 * time: 0.5403568744659424 14 8.290599e+02 6.668149e+00 * time: 0.5404748916625977 15 8.289913e+02 6.310546e+00 * time: 0.5406019687652588 16 8.289476e+02 3.243392e+00 * time: 0.5407168865203857 17 8.289357e+02 6.304275e-01 * time: 0.5408289432525635 18 8.289333e+02 4.855263e-01 * time: 0.5409388542175293 19 8.289316e+02 6.221412e-01 * time: 0.5410749912261963 20 8.289261e+02 1.354835e+00 * time: 0.5411930084228516 21 8.289128e+02 2.424454e+00 * time: 0.5413088798522949 22 8.288769e+02 4.217283e+00 * time: 0.5414290428161621 23 8.287853e+02 7.038105e+00 * time: 0.5415558815002441 24 8.285541e+02 1.137840e+01 * time: 0.5416760444641113 25 8.280335e+02 1.703129e+01 * time: 0.5418078899383545 26 8.272601e+02 2.001094e+01 * time: 0.5419330596923828 27 8.266878e+02 1.255545e+01 * time: 0.5420570373535156 28 8.263927e+02 5.142171e-01 * time: 0.542180061340332 29 8.263814e+02 8.382867e-02 * time: 0.5423159599304199 30 8.263812e+02 2.717007e-02 * time: 0.5424349308013916 31 8.263812e+02 8.128085e-04 * time: 0.5425560474395752
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 64
Observation records: Active Missing
dv: 235 0
Total: 235 0
Number of parameters: Constant Optimized
0 4
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -826.38119
------------------
Estimate
------------------
tvE0 14.485
tvEmax 40.351
tvEC50 207.26
σ_add 8.1464
------------------emax_coefs = coef(emax_fit)
conc_range_ast = range(0, maximum(ast_df.conc), length = 200)
pred_ast = emax_coefs.tvE0 .+ emax_coefs.tvEmax .* conc_range_ast ./ (emax_coefs.tvEC50 .+ conc_range_ast)
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Concentration (ng/mL)", ylabel = "AST")
scatter!(ax, ast_df.conc, ast_df.ast, color = :steelblue, markersize = 5, alpha = 0.5)
lines!(ax, collect(conc_range_ast), pred_ast, color = :red, linewidth = 2, label = "Emax Fit")
axislegend(ax, position = :rb)
fig
Dependent Variable Types:
Data Terminology:
| Method | Description |
|---|---|
| Naive Pooled | Treat all data as from single individual |
| Naive Averaged | Average concentrations across individuals at each time |
| Standard Two-Stage | Individual fits, then summarize across subjects |
| MAP Bayesian | Empirical Bayes estimates with population prior |
| Nonlinear Mixed-Effects | Simultaneous estimation of fixed + random effects |
A typical population PK model has two levels of random effects:
Level 1 — Inter-Individual Variability: \[CL_i = TVCL \cdot e^{\eta_{CL,i}}\]
Level 2 — Residual Variability: \[y_{ij} = f_{ij}(\theta, \eta_i) + \varepsilon_{ij}\]
where \(\eta \sim N(0, \omega^2)\) and \(\varepsilon \sim N(0, \sigma^2)\).
The random effect \(\eta_i\) captures the deviation of individual \(i\) from the population mean. By construction, \(\eta \sim \mathcal{N}(0, \omega^2)\), so its mean is zero — it represents a deviation, not an absolute value. A critical principle: \(\eta\) will be whatever it needs to be to match the observed data. The optimizer does not have knowledge of the underlying physiology; it simply assigns each subject the \(\eta\) that minimizes the objective function, even if that value falls far outside the assumed normal distribution.
The three common IIV parameterizations differ in what they imply about the distribution of the parameter itself:
| Parameterization | Equation | Implied Distribution of \(CL_i\) |
|---|---|---|
| Additive | \(CL_i = TVCL + \eta_i\) | Normal — allows \(CL_i < 0\) |
| Proportional | \(CL_i = TVCL \cdot (1 + \eta_i)\) | Normal — allows \(CL_i < 0\) for large negative \(\eta\) |
| Exponential | \(CL_i = TVCL \cdot e^{\eta_i}\) | Log-normal — constrains \(CL_i > 0\) |
The exponential parameterization is preferred in practice because physiological parameters such as clearance and volume are strictly positive. The log-normal distribution enforced by the exponential form guarantees this constraint without requiring additional boundary handling in the optimizer.
| Model | Equation | Distribution |
|---|---|---|
| Additive | \(CL_i = TVCL + \eta_i\) | Normal CL |
| Proportional | \(CL_i = TVCL \cdot (1 + \eta_i)\) | Normal CL |
| Exponential | \(CL_i = TVCL \cdot e^{\eta_i}\) | Log-Normal CL |
The exponential IIV model (\(CL_i = TVCL \cdot e^{\eta_i}\)) is the most common in population PK because it constrains CL > 0. There are five mathematically equivalent ways to code this in Pumas:
CL = tvCL * exp(η) where η ~ Normal(0, ω²)CL = tvCL * exp(η) where η ~ LogNormal(0, ω²) — Pumas handles the transformCL = exp(log(tvCL) + η) where η ~ Normal(0, ω²)logCL = tvlogCL + η, then CL = exp(logCL) — useful when literature reports log-scale parametersCL ~ LogNormal(log(tvCL), ω) in @randomAll five produce identical distributions. Choose based on:
| Model | Equation | Pumas @derived |
|---|---|---|
| Additive | \(Y_{ij} = F_{ij} + \varepsilon\) | dv ~ Normal(pred, sigma_add) |
| Proportional | \(Y_{ij} = F_{ij} \cdot (1 + \varepsilon)\) | dv ~ Normal(pred, abs(pred) * sigma_prop) |
| Combined | \(Y_{ij} = F_{ij}(1+\varepsilon_1) + \varepsilon_2\) | dv ~ Normal(pred, sqrt(sigma_add^2 + (pred*sigma_prop)^2)) |
| NONMEM | Pumas | Description |
|---|---|---|
| METHOD=0 (FO) | Pumas.FO() |
First-Order, \(\eta=0\) |
| METHOD=1 (FOCE) | Pumas.FOCE() |
First-Order Conditional, \(\eta=\hat\eta_i\) |
| METHOD=1 INTER | Pumas.FOCEI() |
FOCE with \(\eta\)-\(\varepsilon\) interaction |
| METHOD=1 LAPLACIAN | Pumas.LaplaceI() |
Laplacian approximation |
In Pumas, FOCEI is the recommended default for most population models. Use LaplaceI() for models with non-normal likelihoods (e.g., discrete, time-to-event).
Beyond the classical NONMEM-equivalent methods, Pumas provides:
| Method | Pumas | When to Use |
|---|---|---|
| MAP Bayesian | Pumas.MAP() |
Individual estimation with informative priors |
| Full Bayesian MCMC | Pumas.BayesMCMC() |
Small datasets, prior information, uncertainty quantification |
| Variational EM (VEM) | Pumas.VEM() |
Discrete/censored/TTE likelihoods, large populations |
VEM uses variational inference to approximate the posterior — it avoids the inner/outer loop of FOCE and handles non-standard likelihoods natively. See the VEM tutorial for details.
Key diagnostics include:
infer())The complete Pumas workflow follows a systematic pipeline:
@model macro with @param, @random, @pre, @dynamics, @derivedread_pumas() to create a Population from a DataFramefit(model, population, init_params, algorithm) → returns fitted modelinfer(fit) → confidence intervals, RSE, standard errors via coeftable()inspect(fit) → residuals, predictions, individual estimatesgoodness_of_fit(inspect_result) → standard 4-panel GOF plotmetrics_table(fit1, fit2) → AIC, BIC, -2LL side by side; lrtest(fit1, fit2) for nested modelspredict() / simobs() → VPCs via vpc() + vpc_plot()Each step builds on the previous. The inspect() result is reusable across diagnostics, GOF plots, and individual fits — compute it once and pass it around.
In Pumas, inspect() automatically computes the appropriate weighted residuals based on the estimation method used.
# Compare additive vs proportional error for dQTc linear model
# Proportional error model
linear_dqtc_prop = @model begin
@param begin
tvINT ∈ RealDomain(init = 1.0)
tvSLP ∈ RealDomain(init = 0.02)
σprop ∈ RealDomain(lower = 0.0, init = 0.3)
end
@covariates conc
@pre begin
INT = tvINT
SLP = tvSLP
end
@derived begin
mu := @. INT + SLP * conc
dv ~ @. Normal(mu, abs(mu) * σprop + 1e-6)
end
end
prop_lin_params = (tvINT = 1.0, tvSLP = 0.02, σprop = 0.3)
prop_lin_fit = fit(linear_dqtc_prop, dqtc_pop, prop_lin_params, Pumas.NaivePooled())
println("AIC (Additive): ", aic(lin_fit))
println("AIC (Proportional): ", aic(prop_lin_fit))
println("\nLower AIC = better fit to the data")[ 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 8.779226e+03 1.615111e+04 * time: 2.09808349609375e-5 1 3.470889e+03 6.922687e+02 * time: 0.3934648036956787 2 3.416376e+03 7.312349e+02 * time: 0.3936958312988281 3 3.299995e+03 8.429077e+02 * time: 0.3938779830932617 4 3.174559e+03 1.012577e+03 * time: 0.39410996437072754 5 3.034735e+03 1.278121e+03 * time: 0.39432597160339355 6 2.833749e+03 1.860384e+03 * time: 0.39453792572021484 7 2.692216e+03 2.477865e+03 * time: 0.39477992057800293 8 2.453632e+03 4.178792e+03 * time: 0.39502692222595215 9 2.025915e+03 1.195211e+04 * time: 0.39530491828918457 10 1.930324e+03 1.341220e+04 * time: 0.39569997787475586 11 1.883905e+03 5.396638e+03 * time: 0.3960888385772705 12 1.880606e+03 1.088440e+03 * time: 0.3962578773498535 13 1.872726e+03 2.879915e+03 * time: 0.39640283584594727 14 1.834107e+03 1.689536e+04 * time: 0.39654994010925293 15 1.817574e+03 6.910219e+03 * time: 0.3966848850250244 16 1.812441e+03 1.270558e+03 * time: 0.39681482315063477 17 1.807117e+03 2.445604e+03 * time: 0.39694690704345703 18 1.794540e+03 7.405662e+03 * time: 0.3970818519592285 19 1.780023e+03 1.066841e+04 * time: 0.3972148895263672 20 1.765914e+03 1.244111e+04 * time: 0.39734601974487305 21 1.744088e+03 1.286119e+04 * time: 0.3974759578704834 22 1.726898e+03 8.219097e+03 * time: 0.3976438045501709 23 1.724030e+03 5.792934e+03 * time: 0.3978109359741211 24 1.721991e+03 1.579191e+03 * time: 0.3979790210723877 25 1.721060e+03 1.830487e+03 * time: 0.398162841796875 26 1.720636e+03 2.252779e+03 * time: 0.39832091331481934 27 1.720071e+03 1.384128e+01 * time: 0.39846181869506836 28 1.720033e+03 2.025641e+01 * time: 0.3985908031463623 29 1.720031e+03 4.867382e+01 * time: 0.39872288703918457 30 1.720031e+03 5.000592e+00 * time: 0.39885401725769043 31 1.720031e+03 1.181640e-01 * time: 0.3989880084991455 32 1.720031e+03 4.234472e-04 * time: 0.39911794662475586 AIC (Additive): 3572.4847816578085 AIC (Proportional): 3446.0617056462706 Lower AIC = better fit to the data
The following tutorials from PumasTutorials.jl expand on topics introduced in this chapter:
| Topic | Tutorial | Key Additions |
|---|---|---|
| Lecture: PK Terminology & Data | PK W1a | Variables, parameters, constants, and the concept of error in PK data |
| Lecture: Modeling Philosophy | PKPD W1 | Two pillars of pharmacometrics, parsimony, asking the right question |
| Lecture: Structural Models & Estimation | PKPD W2 | Model building principles, estimation approaches, FO/FOCE/Laplacian |
| Lecture: PopPK Model Components & NLME | PKPD W4 | Structural/variability/covariate submodels, η and ε distributions |
| Lecture: 1-Cpt IV Bolus Model | PK W2 | Monoexponential decline, ke, half-life, the model behind the code |
| End-to-end PK workflow | Introduction | NCA → 1-cpt → 2-cpt → model comparison → VPC, all in one analysis |
| Pumas model syntax deep dive | Fitting | Domain types (RealDomain, PDiagDomain, PSDDomain), @param/@random/@pre/@dynamics/@derived block details |
| Dynamical system specification | Dynamical Systems | 9 equivalent ways to code the same model: closed-form, explicit ODE, ModelingToolkit, Catalyst reactions — with linearity detection |
| Log-normal vs. normal parameterization | Log vs Normal | 5 equivalent Pumas code forms for exponential IIV, with simulation proof of equivalence |
| NONMEM OFV comparison | Comparing NONMEM & Pumas | Exact relationship: -2*loglikelihood(fit) == OFV_with_constant; adjustment formula and troubleshooting |