Model Introduction

Understanding the HBV QSP Model and Treatment Mechanisms

Authors

Vijay Ivaturi

David Müller-Widmann

NoteHBV Tutorial Series

This is Tutorial 1 of 6 in the HBV series.

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

1 Learning Outcomes

Hepatitis B virus (HBV) infection is a major global health challenge, with over 290 million people chronically infected. Mathematical models help us understand treatment mechanisms and predict clinical outcomes.

In this tutorial, you will learn:

  • The biology of HBV infection and why it’s difficult to cure
  • The structure of the 11-ODE QSP model capturing HBV dynamics
  • The treatment mechanisms of NA (nucleoside analogs) and IFN (interferon)
  • How to define and simulate the HBV model in Pumas
  • Key clinical endpoints: Functional Cure and HBsAg Loss
TipPrerequisites

This tutorial assumes basic familiarity with:

  • Pumas model syntax (@param, @dynamics, @derived)
  • Ordinary differential equations (ODEs)
  • Pharmacokinetic/pharmacodynamic modeling concepts

If you’ve completed the Tumor Burden Tutorial Series, you have all the background needed.

2 Background

2.1 The HBV Infection Challenge

Hepatitis B Virus (HBV) infects liver cells (hepatocytes) and establishes chronic infection in many patients. Unlike hepatitis C, which can now be cured with antivirals, HBV is much harder to eliminate because:

  1. cccDNA persistence: The virus inserts its genetic material (cccDNA) into the liver, creating a reservoir that survives treatment
  2. Immune exhaustion: Chronic infection leads to T-cell exhaustion, reducing the body’s ability to clear infected cells
  3. HBsAg shedding: Infected cells produce large amounts of surface antigen (HBsAg) that may “decoy” the immune system

2.2 Treatment Goals

Current treatments aim for different levels of control:

Endpoint Definition Clinical Meaning
Viral Suppression HBV DNA < 10^1.4 copies/mL Virus under control (but not eliminated)
HBsAg Loss HBsAg < 0.05 IU/mL Surface antigen cleared
Functional Cure HBsAg Loss + Viral Suppression for ≥24 weeks off-treatment Durable control without ongoing treatment
ImportantWhy Functional Cure Matters

Current treatments (NA therapy) can suppress virus indefinitely, but patients must take medication for life. Functional cure means patients can stop treatment without viral rebound— the “holy grail” of HBV therapy.

2.3 Treatment Modalities

Nucleoside/Nucleotide Analogs (NAs):

  • Examples: Entecavir, Tenofovir
  • Mechanism: Block viral replication by inhibiting reverse transcriptase
  • Effect: Rapidly suppress HBV DNA, slowly reduce HBsAg
  • Limitation: Rarely achieve functional cure alone

Interferon (IFN):

  • Examples: Pegylated interferon-alpha (PEG-IFN)
  • Mechanism: Enhances immune response, promotes infected cell killing
  • Effect: Boosts T-cell activity, accelerates HBsAg decline
  • Limitation: Side effects limit duration, variable response

3 Libraries

using Pumas
using DataFramesMeta
using CategoricalArrays
using AlgebraOfGraphics
using CairoMakie
using Random

4 Part 1: Model Structure Overview

4.1 The 11-ODE System

The HBV QSP model (Cortés‐Ríos et al., 2025) tracks 11 state variables organized into three compartments:

flowchart LR
    subgraph Hepatocytes ["Hepatocyte Populations"]
        T["T<br/>Target cells"]
        I["I<br/>Infected cells"]
        R["R<br/>Resistant cells"]
    end

    subgraph Viral ["Viral Dynamics"]
        V["V<br/>Virus (HBV DNA)"]
        S["S<br/>HBsAg"]
    end

    subgraph Immune ["Immune Response"]
        E["E<br/>Effector T-cells"]
        P["P<br/>Dendritic cells"]
        Q["Q<br/>Delayed signal"]
        X["X<br/>Cytotoxic effect"]
        D["D<br/>Dead cell marker"]
        A["A<br/>ALT proxy"]
    end

    %% Mass flow (solid): production, conversion, transfer
    T -->|"infection"| I
    I -->|"converts"| R
    R -->|"reverts<br/>ρ_o_ρ_r"| T
    I -->|"virus production<br/>p_V×ϵ_NA×ϵ_IFN"| V
    I -->|"HBsAg production<br/>p_S"| S
    D -->|"drives"| A

    %% Modulation (dashed): catalysis, activation, inhibition
    V -.->|"drives<br/>infection"| I
    E -.->|"kills"| I
    E -.->|"kills (m/n)"| T
    E -.->|"kills (m/n)"| R
    X -.->|"promotes<br/>conversion"| R
    E -.->|"produces<br/>debris"| D
    I -.->|"activates"| P
    P -.->|"stimulates"| E
    P -.->|"delays into"| Q
    Q -.->|"exhausts"| E
    I -.->|"activates"| X
    S -.->|"inhibits<br/>(via HBsAg)"| X
    V -.->|"inhibits<br/>(via HBsAg)"| X

    style T fill:#ccffcc,stroke:#00cc00
    style I fill:#ffcccc,stroke:#cc0000
    style R fill:#ffffcc,stroke:#cccc00
    style V fill:#cce5ff,stroke:#0066cc
    style S fill:#cce5ff,stroke:#0066cc
    style E fill:#ffccff,stroke:#cc00cc

4.2 State Variables Explained

State Symbol Description Units
Uninfected hepatocytes T Uninfected liver cells cells
Infected hepatocytes I Cells producing virus cells
Refractory hepatocytes R Cells resistant to cytotoxic effect cells
Viral load V HBV DNA copies/mL
Subviral particles containing surface antigen S HBsAg molecules molecules
Effector T-cells E Immune cells that kill infected cells cells
Dendritic cells P Antigen-presenting cells cells
Delayed signal Q Immune activation delay -
Cytotoxic effect X Cumulative killing signal -
Dead cell marker D Products from cell death -
ALT proxy A Liver inflammation marker -

4.3 Parameter Structure

The model has 23 fixed parameters (from literature/prior studies) and 9 estimated parameters (with inter-individual variability):

Estimated Parameters (with IIV):

Parameter Description Typical Value (log10)
p_S HBsAg production rate 7.643
β (beta) HBV infectivity -4.358
m Immune killing rate 1.618
k_P Antigen presenting cell synthesis -3.968
k_A ALT production rate -4.869
ϵ_NA NA efficacy on HBV production -1.485
ϵ_IFN IFN efficacy on HBV production -1.466
r_E_IFN IFN immune boost 0.313
conv_E Effector T-cell scaling 2.745
NoteLog10-Scale Parameterization

Parameters are on the log10 scale with additive random effects. This achieves log-normal behavior when used in dynamics as 10^param:

Individual value = 10^(typical_value + η)

This is equivalent to multiplicative IIV on the natural scale.

5 Part 2: The Pumas Model Definition

Let’s build the HBV model.

5.1 Model

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)
        # In the latter case, the system is forced to reach a steady state in which
        # I = 0, viral load V = 0, HBsAg S = 0, and antigen presenting cells P = 0.
        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(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
end
PumasModel
  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, I_threshold, σ_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, Ieff, Peff
  Observed: log10_HBsAg, log10_HBV, log10_ALT, log10_CD8, Ieff, Peff

5.2 Parameters

The following parameters are fixed and not varied between patients:

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,
)

The estimates of the remaining parameters:

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,
)

Taken together, the full set of model parameters is

hbv_params = merge(hbv_fixed_params, hbv_estimated_params)

6 Part 3: Understanding Treatment Mechanisms

6.1 Nucleoside Analog Effect

Nucleoside analogs (NAs) block viral production by inhibiting reverse transcriptase. During treatment with NAs (NA = true), \(\epsilon_{NA} < 1\) and hence the viral production is reduced:

\[V' = p_V \cdot \underbrace{\epsilon_{NA}}_{\text{treatment reduction}} \epsilon_{IFN} \cdot I - d_V \cdot V\]

6.2 Interferon Effect

Interferon (IFN) has dual mechanisms:

  1. Viral suppression (like NA):
    • Reduces viral production via \(\epsilon_{IFN}\)
  2. Immune enhancement:
    • Boosts effector T-cell killing via r_E_IFN: \(\text{killing rate} = m \cdot \underbrace{r\_E\_IFN}_{\text{IFN boost}} \cdot E \cdot I\)

6.3 Treatment Arms

The model supports four treatment scenarios:

Arm NA IFN
Control (no treatment) false false
NA only true false
IFN only false true
NA + IFN combination true true

7 Part 4: Basic Simulation

Let’s simulate a typical patient to see the natural disease progression. The model starts from a small viral inoculum and establishes chronic infection over 5 years, mimicking the clinical course from initial infection to chronic steady state.

7.1 Create a Subject

Create a subject with no treatment:

subject_no_tx = Subject(
    id = "No Treatment",
    covariates = (NA = false, IFN = false),
)
Subject
  ID: No Treatment
  Covariates: NA, IFN

7.2 Simulate Natural History

sim_no_tx = simobs(
    hbv_model,
    subject_no_tx,
    hbv_params,
    zero_randeffs(hbv_model, subject_no_tx, hbv_params);
    obstimes = 0:7:(5 * 365),
    simulate_error = false
)

sim_df = DataFrame(sim_no_tx)

7.3 Visualize Natural History

# Prepare data in long format for AoG
sim_long = @chain sim_df begin
    @select :time :log10_HBV :log10_HBsAg :log10_ALT :log10_CD8
    stack(_, Not("time"); variable_name = "biomarker")
    @rsubset! isfinite(:value)
end

# Create plot with AoG
specs = data(sim_long) *
    mapping(:time => (t -> t / 365) => "Time [year]", :value => "log₁₀ concentration") *
    mapping(layout = :biomarker => renamer("log10_HBV" => "HBV DNA [copies/mL]", "log10_HBsAg" => "HBsAg [IU/mL]", "log10_ALT" => "ALT [U/L]", "log10_CD8" => "CD8+ T Cells [cells/mL]")) *
    visual(Lines)
draw(specs; facet = (; linkyaxes = false))

HBV natural history without treatment over 5 years. Top: Viral load (HBV DNA) rises to chronic levels. Bottom: HBsAg reaches a chronic set point. Both remain elevated without intervention.

HBV natural history without treatment over 5 years. Top: Viral load (HBV DNA) rises to chronic levels. Bottom: HBsAg reaches a chronic set point. Both remain elevated without intervention.

8 Part 5: Treatment Scenarios

8.1 Compare Treatment Scenarios

The key to meaningful treatment comparison is a multi-phase clinical protocol: the virus must first establish chronic infection before treatment begins.

  1. Untreated phase (5 years): Virus establishes chronic steady state
  2. NA suppression phase (4 years): Treatment with NAs
  3. Treatment phase (48 weeks): Arm-specific intervention
  4. Off-treatment phase (24 weeks): Assess durability of response
t_chronic = 5 * 365  # 5 years: establish chronic infection
t_na_suppressed = t_chronic + 4 * 365 # 4 years: NA treatment
t_eot = t_na_suppressed + 48 * 7 # end of 48-week treatment
t_end = t_eot + 24 * 7 # end of 24-week follow-up
# Create subjects with time-varying covariates for each treatment arm
treatments = [
    ("No treatment", false, false),
    ("NA", true, false),
    ("IFN", false, true),
    ("NA + IFN", true, true),
]
pop = map(treatments) do (id, tx_NA, tx_IFN)
    return Subject(;
        id,
        covariates_time = [0.0, t_chronic, t_na_suppressed, t_eot],
        covariates = (
            NA = [false, true, tx_NA, false],   # treatment only during Phase 2
            IFN = [false, false, tx_IFN, false],
        ),
    )
end

# Simulate all arms for a "typical" patient (η = 0)
sims = simobs(
    hbv_model,
    pop,
    hbv_params,
    zero_randeffs(hbv_model, pop, hbv_params);
    obstimes = t_na_suppressed:t_end,
    simulate_error = false,
)

# Combine into dataframe
df_sims = stack(DataFrame(sims), ["log10_HBV", "log10_HBsAg", "log10_ALT", "log10_CD8"], ["time", "id"]; variable_name = "biomarker")

# Plot simulations
specs_data = data(df_sims) * mapping(
    :time => (t -> (t - t_na_suppressed) / 7) => "Time after treatment start [week]",
    color = :id => sorter("No treatment", "NA", "IFN", "NA + IFN") => "Treatment",
) * visual(Lines; color = (:blue, 0.2))
specs_refs = data(
    [
        (; value = 1.4, biomarker = "log10_HBV"),
        (; value = log10(0.05), biomarker = "log10_HBsAg"),
    ]
) * mapping(
    linestyle = direct("Functional cure") => "Threshold" => AlgebraOfGraphics.scale(:RefLines)
) * visual(HLines)
specs = (specs_data + specs_refs) * mapping(
    :value => "log₁₀ concentration",
    layout = :biomarker => renamer("log10_HBV" => "HBV DNA [copies/mL]", "log10_HBsAg" => "HBsAg [IU/mL]", "log10_ALT" => "ALT [U/L]", "log10_CD8" => "CD8+ T Cells [cells/mL]"),
)
draw(
    specs,
    scales(RefLines = (; palette = [:dash]));
    axis = (; xticks = [0, 24, 48, 72]),
    facet = (; linkxaxes = true, linkyaxes = false),
)

Effect of treatment on biomarkers in chronic NA-suppressed HBV (5 years of untreated chronic infection + 4 years of NA treatment) for a typical patient. NA reduces HBsAg gradually, IFN provides immune-mediated decline, and combination of NA and IFN achieves the steepest reduction. Gray dashed line marks end of treatment (week 48).

Effect of treatment on biomarkers in chronic NA-suppressed HBV (5 years of untreated chronic infection + 4 years of NA treatment) for a typical patient. NA reduces HBsAg gradually, IFN provides immune-mediated decline, and combination of NA and IFN achieves the steepest reduction. Gray dashed line marks end of treatment (week 48).

8.2 Clinical Endpoints

Common clinical endpoints are:

  1. HBsAg loss: HBsAg < 0.05 IU/mL for 24 weeks post-treatment.

  2. Functional cure: Undetectable HBV DNA < 10^1.4 copies/mL and HBsAg < 0.05 IU/mL for 24 weeks post-treatment.

ImportantWhy does the typical patient not achieve functional cure?

Notice that HBsAg barely declines during treatment and rebounds immediately when treatment stops. This is expected behavior — it reflects a key clinical reality of HBV therapy.

HBsAg production is not directly reduced by drugs. Looking at the ODE: \(S' = p_S \cdot I - d_V \cdot S\). Neither NA nor IFN appear in this equation. While HBsAg is always cleared at rate \(d_V\), its production (\(p_S \cdot I\)) depends entirely on the number of infected cells. For HBsAg to decline meaningfully, \(I\) itself must be reduced.

Infected cell clearance requires immune activation. \(I\) decreases through immune killing (\(m \cdot E \cdot I\)), but in chronic infection the effector T cells (\(E\)) are exhausted (high \(Q\) suppresses \(E\)). The typical patient’s immune system cannot clear \(I\) below the critical threshold (\(10^{-4}\)~cells/mL).

Functional cure requires the threshold switch. When \(I\) drops below \(10^{-4}\)~cells/mL, the model transitions to a “cured” state where viral and HBsAg production stops entirely. At population mean parameters, \(I\) never reaches this threshold — which is why the typical patient rebounds after treatment.

Only patients with favorable random effects achieve FC. Patients with stronger immune killing (higher \(m\)), better IFN response (higher \(r^{E}_{IFN}\)), and less exhaustion can push \(I\) below threshold during treatment. This matches clinical observations: FC rates are only ~5–15% with IFN-based regimens. We demonstrate this mechanism below.

9 Part 6: Population Variability Preview

Let’s see how parameter variability affects outcomes using the multi-phase protocol. We use the estimated IIV from the NLME fit to generate virtual patients and simulate 50 patients to capture the range of outcomes.

# Simulate 50 patients with IIV using the estimated Ω
sims = map(1:50) do id
    subj = Subject(;
        id,
        covariates_time = [0.0, t_chronic, t_na_suppressed, t_eot],
        covariates = (
            NA = [false, true, true, false],
            IFN = [false, false, true, false],
        ),
    )

    # Simulate with IIV:
    # Sampling of random effects is repeated until random effects allow for successful simulation of the dynamics
    randeffs = sample_randeffs(hbv_model, subj, hbv_params)
    sim = simobs(
        hbv_model, subj, hbv_params, randeffs;
        obstimes = t_na_suppressed:t_end,
        simulate_error = false
    )
    while !isvalid(sim)
        randeffs = sample_randeffs(hbv_model, subj, hbv_params)
        sim = simobs(
            hbv_model, subj, hbv_params, randeffs;
            obstimes = t_na_suppressed:t_end,
            simulate_error = false
        )
    end

    return sim
end

# Combine into dataframe
df_sims = stack(DataFrame(sims), ["log10_HBV", "log10_HBsAg", "log10_ALT", "log10_CD8"], ["time", "id"]; variable_name = "biomarker")

# Plot simulations
specs_data = data(df_sims) * mapping(
    :time => (t -> (t - t_na_suppressed) / 7) => "Time after treatment start [week]",
    group = :id,
) * visual(Lines; color = (:gray, 0.4))
specs_refs = data(
    [
        (; value = 1.4, biomarker = "log10_HBV"),
        (; value = log10(0.05), biomarker = "log10_HBsAg"),
    ]
) * mapping(
    linestyle = direct("Functional cure") => "Threshold" => AlgebraOfGraphics.scale(:RefLine),
) * visual(HLines)
specs = (specs_data + specs_refs) * mapping(
    :value => "log₁₀ concentration",
    layout = :biomarker => renamer("log10_HBV" => "HBV DNA [copies/mL]", "log10_HBsAg" => "HBsAg [IU/mL]", "log10_ALT" => "ALT [U/L]", "log10_CD8" => "CD8+ T Cells [cells/mL]"),
)

draw(
    specs,
    scales(RefLine = (; palette = [:dash]));
    axis = (; xticks = [0, 24, 48, 72]),
    facet = (; linkxaxes = true, linkyaxes = false),
    legend = (; position = :bottom),
)

Biomarker dynamics with inter-individual variability under NA+IFN treatment. Each line is a virtual patient sampled from the estimated population distribution. Most patients rebound after treatment stops (week 48), but patients with favorable immune parameters show steeper HBsAg decline.

Biomarker dynamics with inter-individual variability under NA+IFN treatment. Each line is a virtual patient sampled from the estimated population distribution. Most patients rebound after treatment stops (week 48), but patients with favorable immune parameters show steeper HBsAg decline.

9.1 Functional Cure Mechanism: Responder vs Non-Responder

To understand how functional cure works in this model, let’s compare two patients with NA+IFN treatment: the typical patient (population mean, η = 0) and a responder with favorable random effects (stronger immune killing and better IFN response).

The key state variable is infected cells (I): when I crosses below the threshold (\(10^{-4}\) cells/mL), the model switches to a “cured” state where viral and HBsAg production ceases permanently.

# Construct a typical patient and a responder
covariates_time = [0.0, t_chronic, t_na_suppressed, t_eot]
covariates = (NA = [false, true, true, false], IFN = [false, false, true, false])
# Typical patient: η = 0
typical_subject = Subject(; id = "Typical", covariates_time, covariates)
typical_randeffs = zero_randeffs(hbv_model, typical_subject, hbv_params)
# Responder: favorable random effects
# For instance, higher immune enhancement
responder_subject = Subject(; id = "Responder", covariates_time, covariates)
responder_randeffs = (; η = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 0.0])

# Simulate both
sims = simobs(
    hbv_model, [typical_subject, responder_subject], hbv_params,
    [typical_randeffs, responder_randeffs];
    obstimes = t_na_suppressed:t_end,
    simulate_error = false,
)

# Combine into DataFrame
df_sims = @chain DataFrame(sims) begin
    @rtransform! :log10_I = log10(:I)
    stack(_, ["log10_I", "log10_HBV", "log10_HBsAg", "log10_ALT", "log10_CD8"], ["time", "id"]; variable_name = "biomarker")
    @rsubset! isfinite(:value)
end

# Plot simulations
specs_data = data(df_sims) * mapping(
    :time => (t -> (t - t_na_suppressed) / 7) => "Time after treatment start [week]",
    color = :id => sorter("Typical", "Responder") => "Patient",
) * visual(Lines)
specs_refs = data(
    [
        (; value = 1.4, biomarker = "log10_HBV"),
        (; value = log10(0.05), biomarker = "log10_HBsAg"),
        (; value = -4, biomarker = "log10_I"),
    ]
) * mapping(
    linestyle = AlgebraOfGraphics.direct(["Functional cure", "Functional cure", "Infected cells"]) => "Threshold" => AlgebraOfGraphics.scale(:RefLines)
) * visual(HLines)
specs = (specs_data + specs_refs) * mapping(
    :value => "log₁₀ concentration",
    layout = :biomarker => renamer("log10_I" => "Infected cells [cells/mL]", "log10_HBV" => "HBV DNA [copies/mL]", "log10_HBsAg" => "HBsAg [IU/mL]", "log10_ALT" => "ALT [U/L]", "log10_CD8" => "CD8+ T Cells [cells/mL]"),
)

draw(
    specs,
    scales(RefLines = (; palette = [:dash, :dot]));
    axis = (; xticks = [0, 24, 48, 72]),
    facet = (; linkxaxes = true, linkyaxes = false),
    legend = (; position = :bottom, nbanks = 2),
)

Functional cure mechanism. The responder’s stronger immune response drives the concentration of infected cells below the threshold, triggering the switch to durable viral control. Thereby, the responder achieves HBsAg loss that persists off-treatment, while the typical patient rebounds.

Functional cure mechanism. The responder’s stronger immune response drives the concentration of infected cells below the threshold, triggering the switch to durable viral control. Thereby, the responder achieves HBsAg loss that persists off-treatment, while the typical patient rebounds.
TipThe Threshold Switch is Key

The responder achieves functional cure because their stronger immune response (higher \(m\) and \(r_{E,IFN}\)) drives infected cells below \(10^{-4}\) cells/mL. At that point, the model’s threshold switch activates:

  • Viral production stops: \(V' = -d_V \cdot V\) (exponential decay)
  • HBsAg production stops: \(S' = -d_V \cdot S\) (exponential decay)
  • Immune activation ceases, preventing overshoot

This creates a stable cured equilibrium — even after treatment stops, the virus cannot re-establish because the immune system maintains control.

In a Virtual Clinical Trial (Tutorial 6), we sample thousands of virtual patients from the estimated population distribution to predict realistic FC rates.

10 Summary

10.1 Key Concepts

  1. HBV infection is difficult to cure due to cccDNA persistence and immune exhaustion

  2. The 11-ODE QSP model captures:

    • Hepatocyte populations (T, I, R)
    • Viral dynamics (V, S)
    • Immune response (E, P, Q, X, D, A)

  3. Treatment mechanisms:

    • NA: Blocks viral production
    • IFN: Enhances immune killing + reduces viral production

  4. Clinical endpoints:

    • HBsAg Loss: HBsAg < 0.05 IU/mL
    • Functional Cure: HBsAg Loss + viral suppression (HBV DNA < 10^1.4 copies/mL)

10.2 Model Complexity Comparison

Aspect Tumor Burden HBV QSP
State variables 2 11
Estimated parameters 3 9
Fixed parameters 0 22
Treatment covariates 1 2
Clinical endpoints Response rate FC, HBsAg Loss

10.3 What’s Next?

In Tutorial 2: Global Sensitivity Analysis, we’ll identify which parameters most influence HBV treatment outcomes using Sobol and eFAST sensitivity analysis.

10.4 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