Extending the Hierarchical Model for Antibody Kinetics

Overview

Incorporates feedback from Dr. Morrison and Dr.Aiemjoy
Focus exclusively on (Teunis and Eijkeren 2016) model
Clarifies model dynamics: growth, clearance, decay
Uses updated parameter notation: $\mu_y$ , $\mu_b$ , $\gamma$ , $\alpha$ , $\rho$
Assumes block-diagonal covariance structure across biomarkers

Within-Host ODE System (Teunis and Eijkeren 2016)

Two-phase within-host antibody kinetics:

$\frac{dy}{dt} = \begin{cases} \mu_y y(t), & t \le t_1 \\ - \alpha y(t)^\rho, & t > t_1 \end{cases} \quad \text{with } \frac{db}{dt} = \mu_b b(t) - \gamma y(t) \qquad(1)$

Initial conditions: $y(0) = y_0$ , $b(0) = b_0$
Key transition: $t_1$ is the time when $b(t_1) = 0$
Derived quantity: $y_1 = y(t_1)$

Closed-Form Solutions

Antibody concentration $y(t)$

$t \le t_1$ :
$y(t) = y_0 e^{\mu_y t}$
$t > t_1$ :
$y(t) = y_1 \left(1 + (\rho - 1)\alpha y_1^{\rho - 1}(t - t_1)\right)^{- \frac{1}{\rho - 1}}$

Pathogen load $b(t)$

$t \le t_1$ :
$b(t) = b_0 e^{\mu_b t} - \frac{\gamma y_0}{\mu_y - \mu_b} \left(e^{\mu_y t} - e^{\mu_b t} \right)$
$t > t_1$ :
$b(t) = 0$

Time of Peak Response

Peak Time $t_1$

$t_1 = \frac{1}{\mu_y - \mu_b} \log \left( 1 + \frac{(\mu_y - \mu_b) b_0}{\gamma y_0} \right) \qquad(2)$

Peak Antibody Level $y_1$

$y_1 = y_0 e^{\mu_y t_1} \qquad(3)$

Parameter Summary

Table 1: Parameter summary for antibody kinetics model.

Symbol	Description
$\mu_y$	Antibody production rate (growth phase)
$\mu_b$	Pathogen replication rate
$\gamma$	Clearance rate (by antibodies)
$\alpha$	Antibody decay rate
$\rho$	Shape of antibody decay (power-law)
$t_1$	Time of peak response
$y_1$	Peak antibody concentration

Note: Only the first 6 are typically estimated. $y_1$ is derived from the ODE solution at $t_1$ .

Model Comparison: (Teunis and Eijkeren 2016) vs This Presentation

Table 2: Comparison of Teunis (2016) model and this presentation’s model assumptions.

Component	(Teunis and Eijkeren 2016)	This Presentation
Pathogen ODE	$\mu_0 b(t) - c y(t)$	$\mu_b b(t) - \gamma y(t)$
Antibody ODE (pre- $t_1$ )	$\mu y(t)$	$\mu_y y(t)$
Antibody ODE (post- $t_1$ )	$- \alpha y(t)^r$	$- \alpha y(t)^\rho$
Antibody growth type	Pathogen-driven	Self-driven exponential
Antibody rate name	$\mu$	$\mu_y$
$t_1$ formula	Uses $\mu_0$ , $\mu$ , $b_0$ , $c$ , $y_0$	Uses $\mu_b$ , $\mu_y$ , etc.

Note:

(Teunis and Eijkeren 2016) uses linear clearance: $c y(t)$ , not bilinear
Antibody production is driven by pathogen $b(t)$
Our model simplifies by assuming self-expanding antibody dynamics

Full Parameter Model (7 Parameters)

Subject-level parameters:

$\theta_{ij} \sim \mathcal{N}(\mu_j,\, \Sigma_j), \quad \theta_{ij} = \begin{bmatrix} y_{0,ij} \\ b_{0,ij} \\ \mu_{b,ij} \\ \mu_{y,ij} \\ \gamma_{ij} \\ \alpha_{ij} \\ \rho_{ij} \end{bmatrix}$

Hyperparameters – Means:

$\mu_j$ : population-level mean vector for biomarker $j$
Prior on $\mu_j$ :

$\mu_j \sim \mathcal{N}(\mu_{\text{hyp},j},\, \Omega_{\text{hyp},j})$

Core Parameters Used for Curve Drawing

In this presentation, we focus on 5 key parameters required to draw antibody curves:

$y_0$ : initial antibody level
$t_1$ : time of peak antibody response
$y_1$ : peak antibody level
$\alpha$ : decay rate
$\rho$ : shape of decay

Note: $t_1$ and $y_1$ are derived from the full model - These 5 are sufficient for prediction and plotting

Classifying Model Parameters ((Teunis and Eijkeren 2016) Structure)

Estimated Parameters (7 total):

Core model parameters (5): $\mu_b$ , $\mu_y$ , $\gamma$ , $\alpha$ , $\rho$
Initial conditions (2): $y_0$ , $b_0$

Derived Quantity (not estimated):

$y_1$ : peak antibody level computed as $y(t_1)$

Time of Pathogen Clearance $t_1$

Definition: $t_1$ is the time when the pathogen is cleared, i.e., $b(t_1) = 0$

Analytic expression:

$t_1 = \frac{1}{\mu_y - \mu_b} \log\left(1 + \frac{(\mu_y - \mu_b)b_0}{\gamma y_0}\right)$

Key observations: $t_1$ depends on $\mu_b$ , $\mu_y$ , $b_0$ , $y_0$ , and $y_1 = y(t_1)$ is computed based on this time point

Why It’s a Seven-Parameter Model

Our model estimates 7 parameters:
- 5 biological parameters: $\mu_b$ , $\mu_y$ , $\gamma$ , $\alpha$ , $\rho$
- 2 initial conditions: $y_0$ , $b_0$
But we often refer to an 8th quantity: $y_1$
So why isn’t $y_1$ a parameter?

Answer: $y_1$ is a computed value, not directly estimated.

Why $y_1$ Is Not Fit Directly

$y_1$ is the antibody level at the time the pathogen is cleared:

$y_1 = y(t_1) \quad \text{where } b(t_1) = 0$

y1y_1 is not an “input” — it is computed from:
- $\mu_y$ , $y_0$ , $b_0$ , $\mu_b$ , $\gamma$
- via solution of ODEs to find $t_1$ and compute $y(t_1)$

In other words: $y_1$ is a derived output, not a fit parameter.

How $y_1$ Is Computed

$y_1$ is computed by solving the ODE system:

$\frac{dy}{dt} = \mu_y y(t), \quad \frac{db}{dt} = \mu_b b(t) - \gamma y(t)$

Evaluate $y(t)$ at $t = t_1$ using ODE solution:

$y_1 = y(t_1; \mu_y, y_0, b_0, \mu_b, \gamma)$

Recap: What We Estimate

Seven model parameters (7-parameter model for full dynamics):

$\mu_b$ , $\mu_y$ , $\gamma$ , $\alpha$ , $\rho$ (biological process)
$y_0$ , $b_0$ (initial state)

Derived quantity:

$y_1 = y(t_1)$ — not directly estimated, computed

5-parameter subset for curve visualization:

$y_0$ , $y_1$ , $t_1$ , $\alpha$ , $\rho$

Hierarchical Bayesian Structure

Individual parameters:

$\theta_{ij} = \begin{bmatrix} y_{0,ij} \\ b_{0,ij} \\ \mu_{b,ij} \\ \mu_{y,ij} \\ \gamma_{ij} \\ \alpha_{ij} \\ \rho_{ij} \end{bmatrix} \sim \mathcal{N}(\mu_j,\, \Sigma_j)$

Hyperparameters:

$\mu_j$ : population-level means (per biomarker $j$ )
$\Sigma_j$ : $7 \times 7$ covariance matrix over parameters

Subject-Level Parameters: $\theta_{ij}$

$\theta_{ij} \sim \mathcal{N}(\mu_j,\, \Sigma_j), \quad \theta_{ij} \in \mathbb{R}^7$

Where:

$\boldsymbol{\theta}_{ij} = \begin{bmatrix} y_{0,ij} \\ b_{0,ij} \\ \mu_{b,ij} \\ \mu_{y,ij} \\ \gamma_{ij} \\ \alpha_{ij} \\ \rho_{ij} \end{bmatrix}, \quad \Sigma_j \in \mathbb{R}^{7 \times 7}$

Each subject $i$ has a unique 7-parameter vector per biomarker $j$ , capturing individual-level variation in antibody dynamics.

Hyperparameters: Priors on Population Means

Population-level means:

$\boldsymbol{\mu}_j \sim \mathcal{N}(\boldsymbol{\mu}_{\mathrm{hyp},j}, \Omega_{\mathrm{hyp},j})$

Interpretation:

$\boldsymbol{\mu}_j$ : average parameter vector for biomarker $j$
$\boldsymbol{\mu}_{\mathrm{hyp},j}$ : prior guess (e.g., vector of zeros)
$\Omega_{\mathrm{hyp},j}$ : covariance matrix encoding uncertainty

Example:

$\boldsymbol{\mu}_{\mathrm{hyp},j} = 0, \quad \Omega_{\mathrm{hyp},j} = 100 \cdot I_7$

Hyperparameters: Priors on Covariance

Covariance across parameters:

$\Sigma_j^{-1} \sim \mathcal{W}(\Omega_j, \nu_j)$

$\Sigma_j$ : variability/covariance in subject-level parameters
$\Omega_j$ : prior scale matrix
$\nu_j$ : degrees of freedom

Example:

$\Omega_j = 0.1 \cdot I_7, \quad \nu_j = 8$

Measurement Error and Precision Priors

Observed antibody levels:

$\log(y_{\text{obs},ij}) \sim \mathcal{N}(\log(y_{\text{pred},ij}), \tau_j^{-1})$

Precision prior:

$\tau_j \sim \text{Gamma}(a_j, b_j)$

$\tau_j$ : shared measurement precision for biomarker $j$
Gamma prior allows flexible noise modeling

Matrix Algebra Computation

Let $K = 7$ (parameters), $J$ biomarkers. Then:

$\Theta_i = \begin{bmatrix} \theta_{i1} & \theta_{i2} & \cdots & \theta_{iJ} \end{bmatrix} \in \mathbb{R}^{K \times J}$

Assume:

$\text{vec}(\Theta_i) \sim \mathcal{N}(\text{vec}(M), \Sigma_K \otimes I_J)$

Matrix Algebra – Simplified Structure

Setup: $\Theta_i \in \mathbb{R}^{7 \times J}$

Model:

$\text{vec}(\Theta_i) \sim \mathcal{N}(\text{vec}(M), \Sigma_K \otimes I_J)$

$\Sigma_K$ : 7×7 covariance (same across biomarkers)
$I_J$ : biomarkers assumed uncorrelated
Block-diagonal covariance

Understanding $\text{vec}(\Theta_i)$

Each $\theta_{ij} \in \mathbb{R}^7$ :

$\theta_{ij} = \begin{bmatrix} y_0 \\ b_0 \\ \mu_0 \\ \mu_1 \\ c \\ \alpha \\ r \end{bmatrix}$

Flattening:

$\text{vec}(\Theta_i) \in \mathbb{R}^{7J \times 1}$

Understanding $\text{vec}(M)$

Let $M = [\mu_1\, \mu_2\, \cdots\, \mu_J] \in \mathbb{R}^{7 \times J}$

Example for $J=3$ :

$M = \begin{bmatrix} \mu_{1,1} & \mu_{1,2} & \mu_{1,3} \\ \mu_{2,1} & \mu_{2,2} & \mu_{2,3} \\ \mu_{3,1} & \mu_{3,2} & \mu_{3,3} \\ \mu_{4,1} & \mu_{4,2} & \mu_{4,3} \\ \mu_{5,1} & \mu_{5,2} & \mu_{5,3} \\ \mu_{6,1} & \mu_{6,2} & \mu_{6,3} \\ \mu_{7,1} & \mu_{7,2} & \mu_{7,3} \end{bmatrix}$

Covariance Structure: $\Sigma_K \otimes I_J$

$\text{Cov}(\text{vec}(\Theta_i)) = \Sigma_K \otimes I_J$

$\Sigma_K$ : parameter covariance matrix
$I_J$ : biomarker-wise independence
Kronecker product yields block-diagonal matrix

Example: Kronecker Product with $K=2$ , $J=3$

Let:

$\Sigma_K = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{bmatrix},\quad I_3 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$

Then:

$\Sigma_K \otimes I_3 \in \mathbb{R}^{6 \times 6}$

Expanded Matrix: $\Sigma_K \otimes I_3$

$\Sigma_K \otimes I_3 = \begin{bmatrix} \sigma_{11} & 0 & 0 & \sigma_{12} & 0 & 0 \\ 0 & \sigma_{11} & 0 & 0 & \sigma_{12} & 0 \\ 0 & 0 & \sigma_{11} & 0 & 0 & \sigma_{12} \\ \sigma_{21} & 0 & 0 & \sigma_{22} & 0 & 0 \\ 0 & \sigma_{21} & 0 & 0 & \sigma_{22} & 0 \\ 0 & 0 & \sigma_{21} & 0 & 0 & \sigma_{22} \end{bmatrix}$

Next Steps: Modeling Correlation Across Biomarkers

Current Limitation:

Biomarkers assumed independent: $I_J$

Planned Extension:

Use full covariance $\Sigma_J$ :

$\text{Cov}(\text{vec}(\Theta_i)) = \Sigma_K \otimes \Sigma_J$

Extending to Correlated Biomarkers

Assume $K=3$ , $J=3$

Define:

$\Sigma_K = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{bmatrix},\quad \Sigma_J = \begin{bmatrix} \tau_{11} & \tau_{12} & \tau_{13} \\ \tau_{21} & \tau_{22} & \tau_{23} \\ \tau_{31} & \tau_{32} & \tau_{33} \end{bmatrix}$

Kronecker Product Structure: $\Sigma_K \otimes \Sigma_J$

$\Sigma_K \otimes \Sigma_J = \begin{bmatrix} \sigma_{11}\Sigma_J & \sigma_{12}\Sigma_J & \sigma_{13}\Sigma_J \\ \sigma_{21}\Sigma_J & \sigma_{22}\Sigma_J & \sigma_{23}\Sigma_J \\ \sigma_{31}\Sigma_J & \sigma_{32}\Sigma_J & \sigma_{33}\Sigma_J \end{bmatrix}$

Now biomarkers and parameters can be correlated.

Expanded Form: $\Sigma_K \otimes \Sigma_J$ (3x3)

The $9 \times 9$ matrix contains all combinations $\sigma_{ab}\tau_{cd}$

Not block-diagonal — includes cross-biomarker correlation

Practical To-Do List (for Chapter 2)

Model Implementation:

Define full $\Sigma_J$ and prior: $\Sigma_J^{-1} \sim \mathcal{W}(\Psi, \nu)$
Implement $\Sigma_K \otimes \Sigma_J$ in JAGS

Simulation + Validation:

Simulate individuals with correlated biomarkers
Fit both block-diagonal and full-covariance models
Compare fit: DIC, WAIC, predictive checks

Teunis, Peter F. M., and J. C. H. van Eijkeren. 2016. “Linking the Seroresponse to Infection to Within-Host Heterogeneity in Antibody Production.” Epidemics 16: 33–39. https://doi.org/10.1016/j.epidem.2016.04.001.