Virtual Clinical Trial Simulation

Multi-Arm HBV Treatment Comparison with Functional Cure Endpoints

Authors

Vijay Ivaturi

David Müller-Widmann

NoteHBV Tutorial Series

This is Tutorial 6 of 6 in the HBV series.

# Tutorial Status
1 Model Introduction ✓ Complete
2 Global Sensitivity Analysis ✓ Complete
3 Structural Identifiability ✓ Complete
4 Copula Vpop Generation ✓ Complete
5 MILP Calibration ✓ Complete
6 VCT Simulation (current) ← You are here

1 Learning Outcomes

This tutorial brings together all concepts from the HBV series to run a comprehensive Virtual Clinical Trial comparing multiple treatment strategies.

In this tutorial, you will learn:

  • How to structure HBV treatment phases (untreated → NA background → treatment → off-treatment)
  • How to simulate multiple treatment arms (Control, NA, IFN, Combination)
  • How to calculate functional cure and HBsAg loss rates
  • How to analyze and visualize VCT results
TipPrerequisites

This tutorial uses concepts and code from all previous HBV tutorials:

2 Libraries

using Pumas
using DataFramesMeta
using HypothesisTests
using StatsBase
using CategoricalArrays
using Copulas
using JuMP
using HiGHS
using AlgebraOfGraphics
using CairoMakie
using SummaryTables

using Random
using Statistics

# Set seed for reproducibility
Random.seed!(1234)

3 Part 1: HBV VCT Structure

3.1 Treatment Phases

A complete HBV VCT simulates multiple phases:

hbv_phases = (
    untreated = 5 * 365,       # 5 years - establish chronic infection
    na_background = 4 * 365,   # 4 years - NA background (suppressed patients)
    treatment = 48 * 7,        # 48 weeks - active treatment
    off_treatment = 24 * 7,    # 24 weeks - durability assessment
)

3.2 Treatment Arms

treatment_arms = DataFrame(
    arm = ["Control", "NA", "IFN", "NA + IFN"],
    NA = [false, true, false, true],
    IFN = [false, false, true, true],
    description = [
        "No treatment during treatment phase",
        "Nucleoside analog only",
        "Interferon only",
        "Combination therapy",
    ]
)

simple_table(treatment_arms, [:arm => "Arm", :NA, :IFN, :description => "Description"])
Arm NA IFN Description
Control false false No treatment during treatment phase
NA true false Nucleoside analog only
IFN false true Interferon only
NA + IFN true true Combination therapy

4 Part 2: Model, Parameters, and Virtual Population

We use the HBV model and parameters from HBV-01 and HBV-05. See those tutorials for details on the model structure and parameter estimation.

4.1 Model

HBV model definition (same as HBV-01 and HBV-05) 
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()
        "Threshold of infected hepatocytes to produce dynamical switch from progressive disease to long-lasting control of HBV infection [cells/mL]"
        I_threshold ∈ RealDomain(; lower = 0.0)

        # 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 + Ieff + 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 * Ieff - ρ_o_ρ_r * R + r_T * R * (1 - (T + Ieff + R) / T_max) -
            d_TI * R - (m / n) * r_E_IFN * E * R

        # Dead cell marker
        D' = m * r_E_IFN * E * Ieff + (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 * Ieff - d_V * S

        # Viral load (NA and IFN reduce production)
        V' = p_V * Ieff * ϵ_NA * ϵ_IFN - d_V * V

        # Dendritic cell activation
        P' = k_P * (Ieff / (ϕ_E + Ieff)) - d_P * P

        # Effector T cells
        E' = (d_E + r_E * E) * Peff - d_E_Q * E * (Q^4 / (ϕ_Q^4 + Q^4)) - d_E * E

        # Delayed effector signal
        Q' = d_Q * (Peff - Q)

        # Cytotoxic effect
        X' = r_X * (1 - X) * (Ieff / (ϕ_E + Ieff)) *
            (1 - S_max * (HBsAg_ug_mL^nS / (ϕ_S^nS + HBsAg_ug_mL^nS))) - d_X * X
    end

    @vars begin
        # Switch between progressive infection (I ≥ I_threshold) and clearance of infection (I < I_threshold)
        Ieff = ifelse(I < I_threshold, zero(I), I)
        Peff = ifelse(I < I_threshold, zero(P), P)

        # 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(HBsAg_ug_mL / 0.98e-3), σ_add_log10_HBsAg, σ_prop_log10_HBsAg)

        "Viral load [copies/mL]"
        log10_HBV ~ @. CombinedNormal(log10(V), σ_add_log10_HBV, σ_prop_log10_HBV)

        "ALT [U/L]"
        log10_ALT ~ @. ProportionalNormal(log10(A), σ_prop_log10_ALT)

        "CD8+ T cell counts [cells/mL]"
        log10_CD8 ~ @. ProportionalNormal(log10(conv_E * E), σ_prop_log10_CD8)
    end
end

4.2 Parameters

Fixed and estimated parameters 
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,
    I_threshold = 1.0e-4,
)

# IIV standard deviations (on logit/log10 scale)
ω_p_S = 0.632
ω_β = 2.489
ω_m = 0.827
ω_k_P = 1.292
ω_k_A = 1.156
ω_ϵ_NA = 1.583
ω_ϵ_IFN = 1.626
ω_r_E_IFN = 2.739
ω_conv_E = 1.248

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([ω_p_S^2, ω_β^2, ω_m^2, ω_k_P^2, ω_k_A^2, ω_ϵ_NA^2, ω_ϵ_IFN^2, ω_r_E_IFN^2, ω_conv_E^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)

4.3 Copula Sampling and Baseline Simulation

We sample correlated random effects using the Gaussian copula (see HBV-04 and HBV-05).

# Spearman correlation matrix (9×9) from @cortes2025HBV Figure S5B
Rho = [
    # p_S     β      m     k_P    k_A    ϵ_NA   ϵ_IFN  r_E_IFN conv_E
    1.0    0.43   0.46   0.51  -0.7   0.58  -0.14  -0.38    0.0;  # p_S
    0.43   1.0    0.3  -0.01  -0.5   0.38  -0.52  -0.21    0.0;  # β
    0.46   0.3   1.0   -0.16  -0.5   0.29   0.07  -0.33    0.0;  # m
    0.51  -0.01  -0.16   1.0   -0.45   0.25   0.14  -0.1    0.0;  # k_P
    -0.7  -0.5  -0.5  -0.45   1.0   -0.41   0.01   0.05    0.0;  # k_A
    0.58   0.38   0.29   0.25  -0.41   1.0   -0.61  -0.36    0.0;  # ϵ_NA
    -0.14  -0.52   0.07   0.14   0.01  -0.61   1.0    0.26    0.0;  # ϵ_IFN
    -0.38  -0.21  -0.33  -0.1   0.05  -0.36   0.26   1.0     0.0;  # r_E_IFN
    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0      1.0;  # conv_E
]

copula = GaussianCopula(Rho)
correlated_randeffs_dist = SklarDist(
    copula, (
        Normal(0, ω_p_S),
        Normal(0, ω_β),
        Normal(0, ω_m),
        Normal(0, ω_k_P),
        Normal(0, ω_k_A),
        Normal(0, ω_ϵ_NA),
        Normal(0, ω_ϵ_IFN),
        Normal(0, ω_r_E_IFN),
        Normal(0, ω_conv_E),
    )
)

The Everest trial enrolls NA-suppressed patients, so baseline HBsAg is the level after 5 years of untreated chronic infection followed by 4 years of NA suppression (Cortés‐Ríos et al., 2025). To reduce computation time in this tutorial, we only generate 500 virtual patients. We resample random effects that lead to integration errors:

# Sample `nvps` virtual patients
nvps = 500
baseline_pop = map(1:nvps) do id
    Subject(;
        id,
        covariates = (; NA = [false, true], IFN = [false, false]),
        covariates_time = [0, hbv_phases.untreated],
    )
end
baseline_vrandeffs = map(baseline_pop) do _
    return (; η = rand(correlated_randeffs_dist))
end

# Compute baseline HBsAg (after 5 yr chronic infection + 4 year NA suppression)
baseline_obstimes = [hbv_phases.untreated + hbv_phases.na_background]
baseline_sims = simobs(hbv_model, baseline_pop, hbv_params, baseline_vrandeffs; simulate_error = false, obstimes = baseline_obstimes)

while true
    # Check whether all simulations succeeded
    baseline_invalid = findall(!isvalid, baseline_sims)
    isempty(baseline_invalid) && break

    # Resample random effects for failed simulation
    pop = @view(baseline_pop[baseline_invalid])
    vrandeffs = @view(baseline_vrandeffs[baseline_invalid])
    map!(vrandeffs, pop) do _
        return (; η = rand(correlated_randeffs_dist))
    end
    baseline_sims[baseline_invalid] = simobs(hbv_model, pop, hbv_params, vrandeffs; simulate_error = false, obstimes = baseline_obstimes)
end

# Extract vector of baseline HBsAg
baseline_log10_HBsAg = postprocess((gen, _) -> only(gen.log10_HBsAg), baseline_sims)

4.4 ILP Calibration to Everest Trial

We calibrate the virtual population to match the baseline HBsAg distribution from the Everest trial (see HBV-05 for the full ILP theory and derivation):

# Everest trial baseline HBsAg distribution
everest_target = DataFrame(
    hbsag_range = ["<100", "[100, 200)", "[200, 500)", "[500, 1000)", "≥1000"],
    log10_lower = [-Inf, log10(100), log10(200), log10(500), log10(1000)],
    log10_upper = [log10(100), log10(200), log10(500), log10(1000), Inf],
    pct = [37.5, 10.9, 19.1, 21.5, 11.0]
)

simple_table(everest_target, [:hbsag_range => "HBsAg", :pct => "Target (%)"])
HBsAg Target (%)
<100 37.5
[100, 200) 10.9
[200, 500) 19.1
[500, 1000) 21.5
≥1000 11 
# Calibration tolerance
epsilon = 0.05

# Create optimization model
ilp = Model(HiGHS.Optimizer)
set_silent(ilp)
set_attribute(ilp, "presolve", "on")
set_time_limit_sec(ilp, 60.0)

# Decision variables
@variable(ilp, x[1:length(baseline_log10_HBsAg)], Bin)

# Objective: maximize number of selected VPs
ntotal = sum(x)
@objective(ilp, Max, ntotal)

# Constraints: distribution matching for each category
for row in eachrow(everest_target)
    pc = row.pct / 100
    nc = sum(x[row.log10_lower .<= baseline_log10_HBsAg .< row.log10_upper])
    @constraint(ilp, nc <= (pc + epsilon) * ntotal)
    @constraint(ilp, nc >= (pc - epsilon) * ntotal)
end
optimize!(ilp)
assert_is_solved_and_feasible(ilp)
solution_summary(ilp)
solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : 2.34000e+02
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 2.34000e+02
│ ├ dual_objective_value : NaN
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 5.07442e-03
  ├ simplex_iterations : 0
  ├ barrier_iterations : -1
  └ node_count         : 1
selected_vps = findall(xi -> value(xi) == 1, x)
calibrated_log10_HBsAg = baseline_log10_HBsAg[selected_vps]
calibrated_vrandeffs = baseline_vrandeffs[selected_vps]

We compare the baseline HBsAg of the virtual population before and after calibration:

df = vcat(
    DataFrame(; log10_HBsAg = baseline_log10_HBsAg),
    DataFrame(; log10_HBsAg = calibrated_log10_HBsAg);
    source = "stage" => ["Original", "Calibrated"],
)
data_specs = data(df) * mapping(
    :log10_HBsAg,
    col = :stage => sorter("Original", "Calibrated"),
) * histogram(; bins = vcat(everest_target.log10_lower, everest_target.log10_upper[end]), normalization = :probability)
target_specs = mapping(
    collect(Iterators.flatten((row.log10_lower, row.log10_upper) for row in eachrow(everest_target))),
    repeat(everest_target.pct; inner = 2) ./ 100,
    group = repeat(1:nrow(everest_target); inner = 2),
) * visual(Lines; label = "Everest target")
draw(
    data_specs + target_specs;
    axis = (; xlabel = "Baseline HBsAg (log₁₀ IU/mL)", ylabel = "Proportion"),
)

Before and after ILP calibration to Everest trial baseline HBsAg distribution. The calibrated population (right) closely matches the target.

Before and after ILP calibration to Everest trial baseline HBsAg distribution. The calibrated population (right) closely matches the target.

5 Part 3: VCT Simulation

We simulate the full VCT for each treatment arm using the Pumas model with time-varying NA and IFN covariates. Each patient goes through the four phases defined in Part 1.

5.1 Run VCT for All Arms

Biomarkers are measured weekly during treatment and post-treatment, and monthly otherwise:

# Measurement times
vct_obstimes = unique!(
    vcat(
        # Monthly observations during untreated + NA background
        0:30:hbv_phases.untreated,
        hbv_phases.untreated .+ (0:30:hbv_phases.na_background),
        # Weekly observations during treatment and post-treatment
        (hbv_phases.untreated + hbv_phases.na_background) .+ (0:7:hbv_phases.treatment),
        (hbv_phases.untreated + hbv_phases.na_background + hbv_phases.treatment) .+ (0:7:hbv_phases.off_treatment),
    )
)

# Covariate time points
vct_covariates_time = cumsum([0, hbv_phases.untreated, hbv_phases.na_background, hbv_phases.treatment])

# Simulate all treatment arms
vct_all_sims = mapreduce(vcat, eachrow(treatment_arms)) do arm
    pop = map(selected_vps) do id
        return Subject(;
            id,
            covariates = (;
                NA = [false, true, arm.NA, false],
                IFN = [false, false, arm.IFN, false],
            ),
            covariates_time = vct_covariates_time,
        )
    end
    sims = simobs(hbv_model, pop, hbv_params, calibrated_vrandeffs; simulate_error = false, obstimes = vct_obstimes)
    filter!(isvalid, sims)

    sims_df = DataFrame(sims)
    sims_df.arm .= arm.arm

    return sims_df
end

6 Part 4: Endpoint Analysis

6.1 Calculate Functional Cure Rates

Functional cure is defined as HBsAg < 0.05 IU/mL and HBV DNA < 10^1.4 copies/mL for 24 weeks post-treatment.

# End of treatment
t_eot = hbv_phases.untreated + hbv_phases.na_background + hbv_phases.treatment
# End of follow-up
t_fu = t_eot + hbv_phases.off_treatment

fc_df = @by vct_all_sims [:arm, :id] @astable begin
    t_eot_fu = findall(t -> t_eot <= t <= t_fu, :time)
    :functional_cure = all(x -> x < log10(0.05), :log10_HBsAg[t_eot_fu]) && all(x -> x < 1.4, :log10_HBV[t_eot_fu])
end

fc_rate_df = @by fc_df :arm :fc_rate = 100 * mean(:functional_cure)
simple_table(fc_rate_df, [:arm => "Arm", :fc_rate => "Functional Cure (%)"])
Arm Functional Cure (%)
Control 0
NA 0.427
IFN 11.1
NA + IFN 20.5

6.2 Confidence Intervals

We compute 95% confidence intervals (Clopper-Pearson) for the complete response rate:

# Calculate all functional cure rates with CI
fc_ci_df = @by fc_df :arm @astable begin
    npatients_total = length(:functional_cure)
    npatients_fc = sum(:functional_cure)
    :rate = 100 * npatients_fc / npatients_total
    :ci = round.(100 .* confint(BinomialTest(npatients_fc, npatients_total); level = 0.95); sigdigits = 3)
end

simple_table(fc_ci_df, [:arm => "Arm", :rate => "Functional Cure (%)", :ci => "95% CI"])
Arm Functional Cure (%) 95% CI
Control 0 (0.0, 1.56)
NA 0.427 (0.0108, 2.36)
IFN 11.1 (7.39, 15.9)
NA + IFN 20.5 (15.5, 26.3)

7 Part 5: Visualization

7.1 HBsAg Dynamics by Arm

log10_HBsAg_time_summary = @chain vct_all_sims begin
    @rsubset !ismissing(:log10_HBsAg)
    @by [:time, :arm] @astable begin
        q05, q25, q50, q75, q95 = quantile(:log10_HBsAg, (0.05, 0.25, 0.5, 0.75, 0.95))
        :median = q50
        :q25 = q25
        :q75 = q75
        :q05 = q05
        :q95 = q95
    end
end
sort!(log10_HBsAg_time_summary, [:time, :arm])

simple_table(first(log10_HBsAg_time_summary, 10), [:time, :arm => "Arm", :median => "Median", :q25, :q75, :q05, :q95])
Time Arm Median q25 q75 q05 q95
0 Control -3.46 -4.45 -2.76 -6.3 -1.83
0 IFN -3.46 -4.45 -2.76 -6.3 -1.83
0 NA -3.46 -4.45 -2.76 -6.3 -1.83
0 NA + IFN -3.46 -4.45 -2.76 -6.3 -1.83
30 Control 5.87 4.19 6.72 -1.1 7.49
30 IFN 5.87 4.19 6.72 -1.1 7.49
30 NA 5.87 4.19 6.72 -1.1 7.49
30 NA + IFN 5.87 4.19 6.72 -1.1 7.49
60 Control 5.74 3.79 6.53 1.42 7.4
60 IFN 5.74 3.79 6.53 1.42 7.4 
# Time summary
specs_time_summary = data(log10_HBsAg_time_summary) *
    mapping(:time => (t -> t / 7) => "Time [week]", color = :arm => sorter("Control", "NA", "IFN", "NA + IFN") => "Arm") * (
    # IQR band layer
    mapping(:q25, :q75) * visual(Band; alpha = 0.3) +
        # Median line layer
        mapping(:median) * visual(Lines)
)
# Reference lines
specs_references = mapping([log10(0.05)], linestyle = ["Functional cure"] => "Threshold") * visual(HLines)

draw(specs_time_summary + specs_references; axis = (; ylabel = "HBsAg (log₁₀ IU/mL)"))

HBsAg dynamics by treatment arm. Lines show median; shaded regions show IQR.

HBsAg dynamics by treatment arm. Lines show median; shaded regions show IQR.

7.2 Functional Cure Rate Comparison

bar_layer = data(fc_ci_df) *
    mapping(
    :arm => sorter("Control", "NA", "IFN", "NA + IFN") => "Arm"
) * (
    mapping(
        :rate,
    ) *
        visual(BarPlot) +
        mapping(:ci => first, :ci => last) * visual(Rangebars; color = :black, linewidth = 1.5, whiskerwidth = 8)
)
draw(bar_layer; axis = (; ylabel = "Functional cure rate (%)"))

Functional cure rates by treatment arm with 95% confidence intervals.

Functional cure rates by treatment arm with 95% confidence intervals.

7.3 Waterfall Plot

# Calculate change from baseline for combination treatment
# Baseline time
t_bl = hbv_phases.untreated + hbv_phases.na_background
# End of treatment
t_eot = t_bl + hbv_phases.treatment
vct_combo_change = @chain vct_all_sims begin
    @rsubset :arm == "NA + IFN"
    @by :id @astable begin
        :change = only(:log10_HBsAg[:time .== t_eot]) - only(:log10_HBsAg[:time .== t_bl])
    end
    @orderby -:change
    @transform! :rank = 1:length(:change)
end
leftjoin!(vct_combo_change, @rsubset(fc_df, :arm == "NA + IFN"); on = "id")
@rtransform! vct_combo_change :response = :functional_cure ? "Functional cure" : (:change < 0 ? "Partial response" : "No response")

# Data
specs_data = data(vct_combo_change) *
    mapping(
    :rank => "Patient (sorted by change from baseline)",
    :change => "Change from baseline HBsAg (log₁₀ IU/mL)",
    color = :response => sorter("Functional cure", "Partial response", "No response") => "Response",
) *
    visual(BarPlot, gap = 0)
draw(
    specs_data,
)

Waterfall plot: HBsAg change from baseline to end of treatment for combination arm.

Waterfall plot: HBsAg change from baseline to end of treatment for combination arm.

8 Summary

8.1 Key Results

# Best arm
fc_ci_df[argmax(fc_ci_df.rate), :arm]
"NA + IFN"
# Best rate
maximum(fc_ci_df.rate)
20.512820512820515

8.2 VCT Workflow Complete

Congratulations! You have completed the HBV Tutorial Series. You now have the skills to:

  1. Understand HBV biology and the QSP model structure
  2. Analyze parameter sensitivity with GSA
  3. Assess identifiability for clinical measurements
  4. Generate virtual populations with copula sampling
  5. Calibrate to clinical distributions with MILP
  6. Run comprehensive VCTs with multiple treatment arms

8.3 References

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