using Pumas
using GlobalSensitivity
using AlgebraOfGraphics
using CairoMakie
using DataFramesMeta
using SummaryTables
using LinearAlgebra
using Statistics
Global Sensitivity Analysis
Identifying Key Drivers of Treatment Response
NoteHBV Tutorial Series
This is Tutorial 2 of 6 in the HBV series.
| # | Tutorial | Status |
|---|---|---|
| 1 | Model Introduction | ✓ Complete |
| 2 | Global Sensitivity Analysis (current) | ← You are here |
| 3 | Structural Identifiability | |
| 4 | Copula Vpop Generation | |
| 5 | MILP Calibration | |
| 6 | VCT Simulation |
1 Learning Outcomes
Global Sensitivity Analysis (GSA) identifies which parameters most influence model outputs, guiding parameter estimation and experimental design (Cortés-Ríos et al., 2025). For the theory behind GSA (Sobol indices, eFAST, variance decomposition), see TB-02: Global Sensitivity Analysis.
In this tutorial, you will learn:
- How to set up and run GSA on the HBV QSP model
- How to use batch simulation for efficient GSA evaluation
- Which parameters drive HBsAg dynamics (Cortés‐Ríos et al., 2025)
- How to interpret sensitivity indices for complex QSP models
2 Libraries
3 Part 1: HBV Model Setup for GSA
3.1 Model
We use the HBV model from HBV-01, with two modifications for GSA: an @observed block for endpoints, and the Ieff/Peff threshold switch removed (matching the MATLAB implementation in Cortés‐Ríos et al. (2025), which uses I and P directly without a discontinuous switch).
hbv_model = @model begin
@param begin
# Estimated parameters (log10 scale)
"Production/synthesis rate of HBsAg subviral particles from infected hepatocytes [molec./cell/d]"
tv_log10_p_S ∈ RealDomain(; lower = 2, upper = 12)
"HBV infectivity [mL/virus/d]"
tv_log10_β ∈ RealDomain(; lower = -10, upper = 0)
"Killing rate of infected hepatocytes by cellular immunity [mL/cells/d]"
tv_log10_m ∈ RealDomain(; lower = -3, upper = 5)
"Synthesis rate of antigen presenting cells [cells/d]"
tv_log10_k_P ∈ RealDomain()
"Synthesis rate of ALT from debris hepatocytes [10^3 U/cells/d]"
tv_log10_k_A ∈ RealDomain()
"PD efficacy of NA on HBV production []"
tv_log10_ϵ_NA ∈ RealDomain(; lower = -5, upper = 0)
"PD efficacy of IFNα on HBV production []"
tv_log10_ϵ_IFN ∈ RealDomain(; lower = -5, upper = 0)
"PD efficacy of IFNα on activation of cellular immunity []"
tv_log10_r_E_IFN ∈ RealDomain(; lower = 0, upper = 4)
"Scaling factor for observed effector CD4/8 T cell numbers []"
tv_log10_conv_E ∈ RealDomain()
# IIV of estimated parameters in log10 scale
Ω ∈ PDiagDomain(9)
# Fixed parameters
"Baseline of HBV [copies/mL]"
ini_V ∈ RealDomain(; lower = 0.0)
"Baseline of ALT [U/L]"
ini_A ∈ RealDomain(; lower = 0.0)
"Max. proliferation rate of uninfected hepatocytes [1/d]"
r_T ∈ RealDomain(; lower = 0.0)
"Max. proliferation rate of cellular immunity [1/d]"
r_E ∈ RealDomain(; lower = 0.0)
"Max. proliferation rate of innate immunity [1/d]"
r_X ∈ RealDomain(; lower = 0.0)
"Fraction of cellular immunity that killing uninfected and refractory cells []"
n ∈ RealDomain(; lower = 0.0)
"Max. rate of conversion of infected cells to refractory cells due to the action of innate immunity [1/d]"
ρ ∈ RealDomain(; lower = 0.0)
"Max. rate of conversion of refractory cells back to target cells in the absence of innate immune responses []"
ρ_o_ρ_r ∈ RealDomain(; lower = 0.0)
"Production/synthesis rate of HBV from infected hepatocytes [virus/cell/d]"
p_V ∈ RealDomain(; lower = 0.0)
"Death rate of infected hepatocytes [1/d]"
d_TI ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of HBV HBsAg subviral particles [1/d]"
d_V ∈ RealDomain(; lower = 0.0)
"Rate of systemic induction/recovery from exhaustion [1/d]"
d_Q ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of innate immunity [1/d]"
d_X ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of cellular immunity [1/d]"
d_E ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of antigen presenting cells [1/d]"
d_P ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of debris [1/d]"
d_D ∈ RealDomain(; lower = 0.0)
"Degradation/clearance of ALT [1/d]"
d_A ∈ RealDomain(; lower = 0.0)
"Max. rate of exhaustion [1/d]"
d_E_Q ∈ RealDomain(; lower = 0.0)
"Hepatocyte carrying capacity in the liver [cells/mL]"
T_max ∈ RealDomain(; lower = 0.0)
"Half-maximal (saturation) parameter for the effect of systemic exhaustion []"
ϕ_Q ∈ RealDomain(; lower = 0.0)
"Half-maximal (saturation) parameter for the activation of cellular immunity [cells/mL]"
ϕ_E ∈ RealDomain(; lower = 0.0)
"Max. inhibitory effect of HBsAg subviral particles on innate immunity []"
S_max ∈ RealDomain(; lower = 0.0)
"Half-maximal (saturation) parameter for inhibitory effect of HBsAg subviral particles on the generation of innate immunity [μg/mL]"
ϕ_S ∈ RealDomain(; lower = 0.0)
"Hill coefficient for describing HBsAg-induced inhibition of innate immunity []"
nS ∈ RealDomain()
# Residual errors
σ_add_log10_HBsAg ∈ RealDomain(lower = 0.0)
σ_prop_log10_HBsAg ∈ RealDomain(lower = 0.0)
σ_add_log10_HBV ∈ RealDomain(lower = 0.0)
σ_prop_log10_HBV ∈ RealDomain(lower = 0.0)
σ_prop_log10_ALT ∈ RealDomain(lower = 0.0)
σ_prop_log10_CD8 ∈ RealDomain(lower = 0.0)
end
@random begin
η ~ MvNormal(Ω)
end
# Time-varying covariates for treatment phases
@covariates NA IFN
@pre begin
# Individual parameters
# Logit-normal distribution of log10(p_S), scaled to range [2, 12]
p_S = exp10(10 * logistic(logit((tv_log10_p_S - 2) / 10) + η[1]) + 2)
# Logit-normal distribution of log10(β), scaled to range [-10, 0]
β = exp10(10 * (logistic(logit(tv_log10_β / 10 + 1) + η[2]) - 1))
# Logit-normal distribution of log10(m), scaled to range [-3, 5]
m = exp10(8 * logistic(logit((tv_log10_m + 3) / 8) + η[3]) - 3)
# Normal distribution of log10(k_P)
k_P = exp10(tv_log10_k_P + η[4])
# Normal distribution of log10(k_A)
k_A = exp10(tv_log10_k_A + η[5])
# Logit-normal distribution of log10(ϵ_NA), scaled to range [-5, 0]
ϵ_NA = exp10(NA * 5 * (logistic(logit(tv_log10_ϵ_NA / 5 + 1) + η[6]) - 1))
# Logit-normal distribution of log10(ϵ_IFN), scaled to range [-5, 0]
ϵ_IFN = exp10(IFN * 5 * (logistic(logit(tv_log10_ϵ_IFN / 5 + 1) + η[7]) - 1))
# Logit-normal distribution of log10(r_E_IFN), scaled to range [0, 4]
r_E_IFN = exp10(IFN * 4 * logistic(logit(tv_log10_r_E_IFN / 4) + η[8]))
# Normal distribution
conv_E = exp10(tv_log10_conv_E + η[9])
end
# Initial conditions
@init begin
V = ini_V
I = d_V * V / p_V
T = T_max - I
S = p_S * I / d_V
A = ini_A
end
@dynamics begin
# Target uninfected hepatocytes
T' = -β * V * T + r_T * T * (1 - (T + I + R) / T_max) -
d_TI * T + ρ_o_ρ_r * R - (m / n) * r_E_IFN * E * T
# Infected hepatocytes
I' = β * V * T + r_T * I * (1 - (T + I + R) / T_max) -
ρ * X * I - d_TI * I - m * r_E_IFN * E * I
# Resistant hepatocytes
R' = ρ * X * I - ρ_o_ρ_r * R + r_T * R * (1 - (T + I + R) / T_max) -
d_TI * R - (m / n) * r_E_IFN * E * R
# Dead cell marker
D' = m * r_E_IFN * E * I + (m / n) * r_E_IFN * E * (T + R) - d_D * D
# Immune response (ALT proxy)
A' = d_A * (ini_A - A) + k_A * D
# HBsAg dynamics
S' = p_S * I - d_V * S
# Viral load (NA and IFN reduce production)
V' = p_V * I * ϵ_NA * ϵ_IFN - d_V * V
# Dendritic cell activation
P' = k_P * (I / (ϕ_E + I)) - d_P * P
# Effector T cells
E' = (d_E + r_E * E) * P - d_E_Q * E * (Q^4 / (ϕ_Q^4 + Q^4)) - d_E * E
# Delayed effector signal
Q' = d_Q * (P - Q)
# Cytotoxic effect
X' = r_X * (1 - X) * (I / (ϕ_E + I)) *
(1 - S_max * (HBsAg_ug_mL^nS / (ϕ_S^nS + HBsAg_ug_mL^nS))) - d_X * X
end
@vars begin
# HBsAg [μg/mL]: Each particle of V and S is assumed to harbor 96 molecules of S envelope protein
f_molec2ug := 24000 * 1.0e6 / 6.022e23
HBsAg_ug_mL := 96 * (V + S) * f_molec2ug
end
# Error models: Values are log10-transformed
@derived begin
"HBsAg [IU/mL]"
log10_HBsAg ~ @. CombinedNormal(log10(max(0, HBsAg_ug_mL) / 0.98e-3), σ_add_log10_HBsAg, σ_prop_log10_HBsAg)
"Viral load [copies/mL]"
log10_HBV ~ @. CombinedNormal(log10(max(0, V)), σ_add_log10_HBV, σ_prop_log10_HBV)
"ALT [U/L]"
log10_ALT ~ @. ProportionalNormal(log10(max(0, A)), σ_prop_log10_ALT)
"CD8+ T cell counts [cells/mL]"
log10_CD8 ~ @. ProportionalNormal(log10(conv_E * max(0, E)), σ_prop_log10_CD8)
end
endPumasModel
Parameters: tv_log10_p_S, tv_log10_β, tv_log10_m, tv_log10_k_P, tv_log10_k_A, tv_log10_ϵ_NA, tv_log10_ϵ_IFN, tv_log10_r_E_IFN, tv_log10_conv_E, Ω, ini_V, ini_A, r_T, r_E, r_X, n, ρ, ρ_o_ρ_r, p_V, d_TI, d_V, d_Q, d_X, d_E, d_P, d_D, d_A, d_E_Q, T_max, ϕ_Q, ϕ_E, S_max, ϕ_S, nS, σ_add_log10_HBsAg, σ_prop_log10_HBsAg, σ_add_log10_HBV, σ_prop_log10_HBV, σ_prop_log10_ALT, σ_prop_log10_CD8
Random effects: η
Covariates: NA, IFN
Dynamical system variables: T, I, R, D, A, S, V, P, E, Q, X
Dynamical system type: Nonlinear ODE
Derived: log10_HBsAg, log10_HBV, log10_ALT, log10_CD8
Observed: log10_HBsAg, log10_HBV, log10_ALT, log10_CD83.2 Parameter Values
hbv_fixed_params = (
ini_V = exp10(-0.481),
ini_A = exp10(1.25),
r_T = 1.0,
r_E = 0.1,
r_X = 1.0,
n = 100,
ρ = 0.001,
ρ_o_ρ_r = 1.0e-5,
p_V = 100.0,
d_TI = 0.0039,
d_V = 0.67,
d_Q = 0.00385,
d_X = 0.2,
d_E = 0.01,
d_P = 0.239,
d_D = 1.66,
d_A = 0.47,
d_E_Q = 0.3,
T_max = 1.36e7,
ϕ_Q = 0.8,
ϕ_E = 100.0,
S_max = 1.0,
ϕ_S = 0.1473,
nS = 0.486,
)hbv_estimated_params = (;
tv_log10_p_S = 7.643,
tv_log10_β = -4.358,
tv_log10_m = 1.618,
tv_log10_k_P = -3.968,
tv_log10_k_A = -4.869,
tv_log10_ϵ_NA = -1.485,
tv_log10_ϵ_IFN = -1.466,
tv_log10_r_E_IFN = 0.313,
tv_log10_conv_E = 2.745,
Ω = Diagonal([0.632^2, 2.489^2, 0.827^2, 1.292^2, 1.156^2, 1.583^2, 1.626^2, 2.739^2, 1.248^2]),
σ_add_log10_HBsAg = 0.0972,
σ_prop_log10_HBsAg = 2.2204e-16,
σ_add_log10_HBV = 0.3206,
σ_prop_log10_HBV = 0.0838,
σ_prop_log10_ALT = 0.0478,
σ_prop_log10_CD8 = 0.4149,
)hbv_params = merge(hbv_fixed_params, hbv_estimated_params)3.3 Base Parameters for GSA
For GSA, we fix random effects and residual error to zero (see TB-02 for why):
base_params = merge(
hbv_params,
(;
Ω = zero(hbv_params.Ω),
σ_add_log10_HBsAg = 0.0,
σ_prop_log10_HBsAg = 0.0,
σ_add_log10_HBV = 0.0,
σ_prop_log10_HBV = 0.0,
σ_prop_log10_ALT = 0.0,
σ_prop_log10_CD8 = 0.0,
)
)3.4 Clinical Protocol
Following Cortés‐Ríos et al. (2025), GSA is performed on NA-suppressed patients — chronic HBV patients who received 4 years of NA background therapy before the treatment period.
t_chronic = 5 * 365 # End of chronic infection phase
t_na_suppressed = t_chronic + 4 * 365 # End of 4-year NA suppression
t_eot = t_na_suppressed + 48 * 7 # End of 48-week treatment
t_end = t_eot + 24 * 7 # End of 24-week off-treatment follow-upWe study the sensitivity of HBsAg levels at different time points: at the initial time point, after 5 years of chronic infection, after 4 years of NA suppression, after 48 weeks of treatment with a combination treatment with IFN and NA, and after 24 weeks of off-treatment follow-up.
gsa_obstimes = [0.0, t_chronic, t_na_suppressed, t_eot, t_end]
gsa_covariates_time = [0.0, t_chronic, t_na_suppressed, t_eot]
gsa_covariates = (; NA = [false, true, true, false], IFN = [false, false, true, false])(NA = Bool[0, 1, 1, 0], IFN = Bool[0, 0, 1, 0])3.5 Parameter Domains
We analyze all 9 estimated parameters. As seen in HBV-01, extreme parameter values can lead to simulation failures. For the global sensitivity analysis, we have to ensure that the model is simulated only at parameter values for which no integration errors occur, as otherwise HBsAg levels are predicted as NaN which contaminates the analysis and yields NaN sensitivity indices.
As a simple heuristic in this tutorial, we choose parameter domains based on the estimated random effect distributions, ensuring that the model can be simulated successfully with both the minimal and maximal values in the parameter domains:
# Lower bound of parameter values
sim_lower = simobs(
hbv_model,
Subject(; time = gsa_obstimes, covariates_time = gsa_covariates_time, covariates = gsa_covariates),
base_params,
(; η = -sqrt.(diag(hbv_params.Ω)) / 3);
simulate_error = false,
)
params_lower = postprocess(sim_lower) do gen, _
idx = only(findall(gen.IFN .&& gen.NA))
return (; tv_log10_p_S = log10(gen.p_S[idx]), tv_log10_β = log10(gen.β[idx]), tv_log10_m = log10(gen.m[idx]), tv_log10_k_P = log10(gen.k_P[idx]), tv_log10_k_A = log10(gen.k_A[idx]), tv_log10_ϵ_NA = log10(gen.ϵ_NA[idx]), tv_log10_ϵ_IFN = log10(gen.ϵ_IFN[idx]), tv_log10_r_E_IFN = log10(gen.r_E_IFN[idx]), tv_log10_conv_E = log10(gen.conv_E[idx]))
end
# Upper bound of parameter values
sim_upper = simobs(
hbv_model,
Subject(; time = gsa_obstimes, covariates_time = gsa_covariates_time, covariates = gsa_covariates),
base_params,
(; η = sqrt.(diag(hbv_params.Ω)) / 3);
simulate_error = false,
)
params_upper = postprocess(sim_upper) do gen, _
idx = only(findall(gen.IFN .&& gen.NA))
return (; tv_log10_p_S = log10(gen.p_S[idx]), tv_log10_β = log10(gen.β[idx]), tv_log10_m = log10(gen.m[idx]), tv_log10_k_P = log10(gen.k_P[idx]), tv_log10_k_A = log10(gen.k_A[idx]), tv_log10_ϵ_NA = log10(gen.ϵ_NA[idx]), tv_log10_ϵ_IFN = log10(gen.ϵ_IFN[idx]), tv_log10_r_E_IFN = log10(gen.r_E_IFN[idx]), tv_log10_conv_E = log10(gen.conv_E[idx]))
endAll other parameters are fixed to their base values:
gsa_constantcoef = (
keys(hbv_fixed_params)...,
:Ω, :σ_add_log10_HBsAg, :σ_prop_log10_HBsAg, :σ_add_log10_HBV, :σ_prop_log10_HBV, :σ_prop_log10_ALT, :σ_prop_log10_CD8,
)(:ini_V, :ini_A, :r_T, :r_E, :r_X, :n, :ρ, :ρ_o_ρ_r, :p_V, :d_TI, :d_V, :d_Q, :d_X, :d_E, :d_P, :d_D, :d_A, :d_E_Q, :T_max, :ϕ_Q, :ϕ_E, :S_max, :ϕ_S, :nS, :Ω, :σ_add_log10_HBsAg, :σ_prop_log10_HBsAg, :σ_add_log10_HBV, :σ_prop_log10_HBV, :σ_prop_log10_ALT, :σ_prop_log10_CD8)3.6 Run eFAST
For demonstration purposes and to keep computation times in this tutorial small, we perform sensitivity analysis with eFAST and with only 75 samples. Moreover, we parallelize simulations using multithreading by specifying batch = true and ensemblealg = EnsembleThreads() keyword arguments:
gsa_efast = gsa(
hbv_model,
Subject(; time = gsa_obstimes, covariates_time = gsa_covariates_time, covariates = gsa_covariates),
base_params,
eFAST(),
[:log10_HBsAg],
params_lower,
params_upper;
samples = 75,
constantcoef = gsa_constantcoef,
ensemblealg = EnsembleThreads(),
batch = true,
)4 Part 3: Interpreting Results
4.1 Extract Indices
We combine first- and total-order indices in a single DataFrame for easier downstream processing:
gsa_efast_first_order = stack(
gsa_efast.first_order,
Not(:dv_name);
variable_name = :parameter,
value_name = :first_order,
)
gsa_efast_total_order = stack(
gsa_efast.total_order,
Not(:dv_name);
variable_name = :parameter,
value_name = :total_order,
)
gsa_efast_df = outerjoin(
gsa_efast_first_order,
gsa_efast_total_order; on = [:dv_name, :parameter],
)
# Compute interaction (difference between total- and first-order indices)
@rtransform! gsa_efast_df :interaction = :total_order - :first_order
# Parse time
@rtransform! gsa_efast_df @astable begin
m = match(r"(?<dv>[^\[]+)\[(?<t>\d+)\]", :dv_name)
:dv = m["dv"]
:time = gsa_obstimes[parse(Int, m["t"])]
end
@select! gsa_efast_df Not(:dv_name)
simple_table(gsa_efast_df, [:time => "Time [d]", :parameter => "Parameter", :first_order => "First order", :total_order => "Total order", :interaction => "Interaction"])4.2 Identify Influential Parameters
threshold = 0.1
influential_gsa_efast_df = @chain gsa_efast_df begin
@rtransform :interpretation = :total_order >= threshold ? "INFLUENTIAL" : "minor"
sort!(_, [:dv, :time, order(:total_order; rev = true)])
end
simple_table(influential_gsa_efast_df, [:dv => "Output", :time => "Time [d]", :parameter => "Parameter", :total_order => "Total order", :interpretation => "Interpretation"])| Output | Time [d] | Parameter | Total order | Interpretation |
| log10_HBsAg | 0 | tv_log10_p_S | 1 | INFLUENTIAL |
| log10_HBsAg | 0 | tv_log10_conv_E | 0.00229 | minor |
| log10_HBsAg | 0 | tv_log10_m | 0.00229 | minor |
| log10_HBsAg | 0 | tv_log10_k_P | 0.00229 | minor |
| log10_HBsAg | 0 | tv_log10_β | 0.00229 | minor |
| log10_HBsAg | 0 | tv_log10_k_A | 0.00229 | minor |
| log10_HBsAg | 0 | tv_log10_ϵ_IFN | 0.00228 | minor |
| log10_HBsAg | 0 | tv_log10_r_E_IFN | 0.00228 | minor |
| log10_HBsAg | 0 | tv_log10_ϵ_NA | 0.00228 | minor |
| log10_HBsAg | 1825 | tv_log10_m | 0.942 | INFLUENTIAL |
| log10_HBsAg | 1825 | tv_log10_k_P | 0.839 | INFLUENTIAL |
| log10_HBsAg | 1825 | tv_log10_p_S | 0.0755 | minor |
| log10_HBsAg | 1825 | tv_log10_β | 0.00398 | minor |
| log10_HBsAg | 1825 | tv_log10_ϵ_NA | 0.00397 | minor |
| log10_HBsAg | 1825 | tv_log10_r_E_IFN | 0.00397 | minor |
| log10_HBsAg | 1825 | tv_log10_ϵ_IFN | 0.00397 | minor |
| log10_HBsAg | 1825 | tv_log10_k_A | 0.00397 | minor |
| log10_HBsAg | 1825 | tv_log10_conv_E | 0.00396 | minor |
| log10_HBsAg | 3285 | tv_log10_k_P | 0.989 | INFLUENTIAL |
| log10_HBsAg | 3285 | tv_log10_m | 0.979 | INFLUENTIAL |
| log10_HBsAg | 3285 | tv_log10_p_S | 0.203 | INFLUENTIAL |
| log10_HBsAg | 3285 | tv_log10_β | 0.0897 | minor |
| log10_HBsAg | 3285 | tv_log10_conv_E | 0.0895 | minor |
| log10_HBsAg | 3285 | tv_log10_k_A | 0.089 | minor |
| log10_HBsAg | 3285 | tv_log10_ϵ_IFN | 0.0886 | minor |
| log10_HBsAg | 3285 | tv_log10_r_E_IFN | 0.0885 | minor |
| log10_HBsAg | 3285 | tv_log10_ϵ_NA | 0.0879 | minor |
| log10_HBsAg | 3621 | tv_log10_β | 0.769 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_m | 0.612 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_k_P | 0.606 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_ϵ_IFN | 0.457 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_r_E_IFN | 0.411 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_ϵ_NA | 0.346 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_p_S | 0.236 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_conv_E | 0.16 | INFLUENTIAL |
| log10_HBsAg | 3621 | tv_log10_k_A | 0.16 | INFLUENTIAL |
| log10_HBsAg | 3789 | tv_log10_m | 0.942 | INFLUENTIAL |
| log10_HBsAg | 3789 | tv_log10_k_P | 0.81 | INFLUENTIAL |
| log10_HBsAg | 3789 | tv_log10_p_S | 0.212 | INFLUENTIAL |
| log10_HBsAg | 3789 | tv_log10_β | 0.196 | INFLUENTIAL |
| log10_HBsAg | 3789 | tv_log10_ϵ_IFN | 0.0535 | minor |
| log10_HBsAg | 3789 | tv_log10_ϵ_NA | 0.0466 | minor |
| log10_HBsAg | 3789 | tv_log10_conv_E | 0.0385 | minor |
| log10_HBsAg | 3789 | tv_log10_k_A | 0.0379 | minor |
| log10_HBsAg | 3789 | tv_log10_r_E_IFN | 0.0105 | minor |
NoteInterpreting HBV GSA Results
Unlike the tumor burden model where a single parameter (f) dominates, the HBV QSP model has a richer sensitivity structure due to its complex interactions between viral dynamics, immune response, and treatment.
Key questions to consider:
- Which parameters are influential (Sₜ > 0.1)? These require careful estimation from clinical data.
- Which parameters show large interactions (Sₜ - S₁ > 0.1)? These suggest nonlinear coupling between model components.
- Which parameters are minor (Sₜ < 0.1)? These can potentially be fixed to literature values.
5 Part 4: Visualization
5.1 Bar Chart: First-Order vs Total-Order
long_df = stack(
gsa_efast_df[gsa_efast_df.time .== gsa_obstimes[end], [:parameter, :first_order, :interaction]],
[:first_order, :interaction];
variable_name = :type,
value_name = :index,
)
parameter_labels = [
"tv_log10_p_S" => rich("p", subscript("S")),
"tv_log10_β" => "β",
"tv_log10_m" => "m",
"tv_log10_k_P" => rich("k", subscript("P")),
"tv_log10_k_A" => rich("k", subscript("A")),
"tv_log10_ϵ_NA" => rich("ϵ", subscript("NA")),
"tv_log10_ϵ_IFN" => rich("ϵ", subscript("IFN")),
"tv_log10_r_E_IFN" => rich("r", subsup("IFN", "E")),
"tv_log10_conv_E" => rich("conv", subscript("E")),
]
specs_gsa = data(long_df) *
mapping(
:parameter => sorter(map(first, parameter_labels)) => "Parameter",
:index => "Sensitivity",
stack = :parameter,
color = :type => renamer(["first_order" => "First-order (Sᵢ)", "interaction" => "Interaction (Sₜᵢ - Sᵢ)"]) => "Type",
) * visual(BarPlot)
specs_threshold = mapping([0.1], color = "Sₜᵢ = 0.1" => "Importance threshold" => AlgebraOfGraphics.scale(:ColorLine)) * visual(HLines; linestyle = :dash)
draw(
specs_gsa + specs_threshold, scales(ColorLine = (; palette = [:red]));
axis = (; xticks = (1:length(parameter_labels), map(last, parameter_labels))), legend = (; position = :bottom)
)
6 Summary
6.1 Key Takeaways
Batch simulation via
simobson a population enables efficient evaluation of the many model runs required for GSAThe HBV QSP model has a richer sensitivity structure than the simpler tumor burden model, with multiple potentially influential parameters
For comprehensive multi-output, multi-timepoint results, see Cortés‐Ríos et al. (2025) Figure S10
For general GSA methodology (sample size selection, interpretation guidelines, common pitfalls), see TB-02.
6.2 Quick Reference: Key Functions
| Function | Package | Purpose |
|---|---|---|
GlobalSensitivity.gsa(f, method, bounds; batch, samples) |
GlobalSensitivity | Run GSA with batch evaluation |
eFAST() |
GlobalSensitivity | eFAST method (fast screening) |
simobs(model, pop, params, vrandeffs) |
Pumas | Parallel population simulation |
6.3 What’s Next?
In Tutorial 3: Structural Identifiability, we’ll analyze which HBV parameters can be identified from different measurement scenarios.
6.4 References
Cortés-Ríos, J., Magee, M., Sher, A., Jusko, W. J., & Desikan, R. (2025). A step-by-step workflow for performing in silico clinical trials with nonlinear mixed effects models. CPT: Pharmacometrics & Systems Pharmacology. https://doi.org/10.1002/psp4.70122
Cortés‐Ríos, J., Ren, T., Hanan, N., Sher, A., Nader, A., Magee, M., Jusko, W. J., & Desikan, R. (2025). Hepatitis b in silico trials capture functional cure, indicate mechanistic pathways, and suggest prognostic biomarker signatures. Clinical Pharmacology &Amp; Therapeutics, 118(3), 600–612. https://doi.org/10.1002/cpt.3718