---
title: "Probabilistic Model of Bilateral Lymphatic Spread in Head and Neck Cancer"
abstract: |
  Purpose: Current guidelines for elective nodal irradiation in oropharyngeal squamous cell carcinoma (OPSCC) recommend including large portions of the contralateral lymphatic system in the clinical target volume (CTV-N), even for lateralized tumors with no clinical lymph node involvement in the contralateral neck. This study introduces a probabilistic model of bilateral lymphatic tumor progression in OPSCC to estimate personalized risks of occult disease in specific lymph node levels (LNLs) based on clinical lymph node involvement, T-stage, and tumor lateralization.

  Methods: Building on a previously developed hidden Markov model for ipsilateral lymphatic spread, we extend the approach to contralateral neck involvement. The model represents LNLs I, II, III, IV, V, and VII on both sides of the neck as binary hidden variables (healthy or involved), connected via arcs representing spread probabilities. These probabilities are learned using Markov chain Monte Carlo (MCMC) sampling from a dataset of 833 OPSCC patients, enabling the model to reflect the underlying lymphatic progression dynamics.

  Results: The model accurately and precisely describes observed patterns of lymph node involvement with a compact set of interpretable parameters. Midline extension of the primary tumor is identified as the primary risk factor for contralateral involvement, with advanced T-stage and extensive ipsilateral involvement further increasing risk. Occult disease in contralateral LNL III is highly unlikely if upstream LNL II is clinically negative, and in contralateral LNL IV, occult disease is exceedingly rare without LNL III involvement.

  Conclusions: This model offers an interpretable, probabilistic framework to inform personalized elective CTV-N volume reduction. For lateralized tumors that do not cross the midline, it suggests the contralateral neck may safely be excluded from elective irradiation. For tumors extending across the midline but with a clinically negative contralateral neck, elective irradiation could be limited to LNL II, reducing unnecessary exposure of normal tissue while maintaining regional tumor control.
---

```{python}
#| echo: false
import re

import pandas as pd
import numpy as np

from lymph import models
from lyscripts import utils
from lyscripts.plot.utils import COLORS

from scripts import shared, paths
from scripts.shared import COL, get_lnl_cols

simple_model = shared.get_model("simple", load_samples=True)
```

```{python}
#| tags: [parameters]
#| echo: false

# width of figures. May depend on number of text columns
full = 17 # cm
half = full
# half =  8 # cm
```

# Introduction

In head and neck squamous cell carcinomas (HNSCC) treatments with radiotherapy or surgery, both the primary tumor and clinically detected lymph node metastases are targeted. In addition, current guidelines include large portions of the neck in the elective clinical target volume (CTV-N) [@gregoire_ctbased_2003;@gregoire_delineation_2014;@gregoire_delineation_2018;@eisbruch_intensitymodulated_2002;@biau_selection_2019;@chao_determination_2002;@vorwerk_guidelines_2011;@ferlito_elective_2009] to mitigate the risk of regional recurrences from untreated microscopic disease undetectable by in-vivo imaging modalities such as computed tomography (CT), magnetic resonance imaging (MRI), or positron emission tomography (PET). However, this approach must balance minimizing the risk of occult disease in the lymphatic drainage region against the toxicity of unnecessarily irradiating healthy tissue.

These CTV-N guidelines rely on anatomically defined lymph node levels (LNLs) [@gregoire_delineation_2014] and the overall prevalence of lymph node metastases within these levels. They often recommend extensive irradiation of both sides of the neck. However, the general prevalence of metastasis in a given LNL does not correspond to an individual patient's risk of occult disease in that region, which depends on their specific state of tumor progression. For example, a patient with no clinically detectable nodal disease (cN0) who has a small, clearly lateralized T1 tumor would receive the same contralateral CTV-N as a patient with significant ipsilateral nodal involvement and an advanced tumor crossing the mid-sagittal plane. Both patients receive elective irradiation of the contralateral LNLs II, III, and IVa [@biau_selection_2019].

To better quantify individualized risk of occult disease, we previously developed an intuitive probabilistic hidden Markov model (HMM) [@ludwig_hidden_2021;@ludwig_modelling_2023], originally based on a conceptually similar a Bayesian network model [@pouymayou_bayesian_2019]. However, these models have been limited to predicting ipsilateral nodal involvement. This work extends the model to include contralateral risk predictions, enabling more personalized radiation volume recommendations for the contralateral neck. By identifying patients with low contralateral risk, the model could guide reductions in the contralateral CTV-N, thereby decreasing radiation-induced toxicity and improving quality of life.

The main contributions of this paper are as follows:

1. @Sec-data presents a multi-centric dataset on lymph node involvement in {{< var data.num_patients >}} OPSCC patients, identifying key risk factors for contralateral lymph node involvement and outlining requirements for a bilateral model extension (@sec-requirements).

2. @Sec-ext-to-contra introduces a bilateral HMM that incorporates primary tumor lateralization, T-category, and clinical involvement as risk factors for contralateral involvement. Model training and computational experiments are described in @sec-methods.

3. @Sec-results demonstrates the model's ability to replicate observed contralateral lymph node involvement patterns and estimates occult disease risk for typical patients. Implications for volume-deescalated radiotherapy are discussed in @sec-discussion.


# Data on Lymphatic Progression Patterns {#sec-data}

To develop models for lymphatic tumor progression for all relevant LNLs, including contralateral regions, we compiled a detailed dataset of {{< var data.num_patients >}} patients with newly diagnosed oropharyngeal squamous cell carcinomas (OPSCC) [@ludwig_dataset_2022;@ludwig_multicentric_2024]. The dataset includes lymph node involvement per LNL for each patient in tabular form, along with primary tumor and patient characteristics such as T-category, subsite, primary tumor lateralization, and HPV p16 status. Patient records were collected from four institutions, and an overview of patient characteristics is provided in @tbl-data-overview.

Data from Inselspital Bern (ISB) and Centre Léon Bérard (CLB) consist exclusively of patients who underwent neck dissections. In contrast, the majority of patients from the University Hospital Zürich (USZ) and the Hospital Vall d'Hebron (HVH) were treated with definitive radiotherapy. Since surgical treatment is more common for early T-category patients, ISB and CLB datasets include a higher proportion of these cases compared to USZ and HVH. For 83 patients in the CLB dataset, the primary tumor's lateralization was not reported.

```{python}
#| echo: false
#| label: tbl-data-overview
#| tbl-cap: Overview over the four datasets from four different institutions used to train and evaluate our model. Here, we briefly characterize the total number of OPSCC patients from the respective institution, their median age, what proportion received neck dissection, the N0 portion of patients, what percentage presented with early T-category (T1/T2), and the prevalence of primary tumor midline extension. For a much more detailed look at the data, visit [lyprox.org](https://lyprox.org).
from typing import Literal

def align(
  column: pd.Series,
  where: Literal["right", "left", "center"],
) -> list[str]:
  """Align a column."""
  return [f"text-align: {where}"] * len(column)

def right_align(column: pd.Series) -> list[str]:
  """Right align column."""
  return align(column, where="right")

def highlight(column: pd.Series, mapping: dict) -> list[str]:
  """Color based on `mapping`."""
  colors = column.map(mapping)
  return colors.map(lambda x: f"color: {x}")

def highlight_bool(
  column: pd.Series,
  c_false: str = COLORS["green"],
  c_true: str = COLORS["red"],
) -> list[str]:
  """Make cell read (`False`) or green (`True`)."""
  return highlight(column, mapping={True: c_true, False: c_false})

def highlight_t_stage(
  column: pd.Series,
  c_early: str = COLORS["green"],
  c_late: str = COLORS["red"],
) -> list[str]:
  """Highlight `early` and `late`."""
  is_early = column == "early"
  return highlight_bool(is_early, c_false=c_late, c_true=c_early)


col_map = {
  COL.inst: "Institution",
  COL.age: "Age",
  COL.nd: "Neck Dissection",
  COL.t_stage: "T-category",
  COL.n_stage: "N-category",
}
short_col_map = {tpl[-1]: val for tpl, val in col_map.items()}

raw = utils.load_patient_data(paths.data)
subdata = raw[col_map.keys()]
subdata.columns = subdata.columns.get_level_values(2)
subdata = (
  subdata
  .rename(columns=short_col_map)
  .reset_index(drop=True)
)

def n0_mean(n_stage: pd.Series) -> float:
  """Compute portion of N0 patients."""
  return (n_stage == 0).mean()

def early_mean(t_stage: pd.Series) -> float:
  """Compute early T-category portion."""
  return t_stage.isin([0,1,2]).mean()

grouped = (
  raw.groupby(
    by=COL.inst,
  ).aggregate(**{
    "Total": pd.NamedAgg(column=COL.age, aggfunc="count"),
    "Age (median)": pd.NamedAgg(column=COL.age, aggfunc="median"),
    "Neck Dissection": pd.NamedAgg(column=COL.nd, aggfunc="mean"),
    "N0": pd.NamedAgg(column=COL.n_stage, aggfunc=n0_mean),
    "Early T-Cat.": pd.NamedAgg(column=COL.t_stage, aggfunc=early_mean),
    "Mid. Ext.": pd.NamedAgg(column=COL.midext, aggfunc="mean"),
  }).convert_dtypes()
)
grouped.index.name = "Institution"
(
  grouped
  .reset_index()
  .style
  .format(
    formatter="{:>.0%}",
    subset=["Neck Dissection", "N0", "Early T-Cat.", "Mid. Ext."],
  )
  .apply(
    func=right_align,
    axis="index",
    subset=grouped.columns[0:],
  )
  .hide()
)
```


## Consensus on Involvement Status {#sec-data-consensus}

Pathological involvement is available only for surgically treated patients and for the levels that were dissected. For non-surgical patients, involvement status is determined clinically, i.e. using imaging. For this work, diagnostic information was synthesized into a consensus decision for each patient and LNL. This consensus reflects the most likely state of involvement and accounts for the sensitivity and specificity of various diagnostic modalities, as reported in the literature [@debondt_detection_2007;@kyzas_18ffluorodeoxyglucose_2008].

The consensus process is detailed in @sec-consensus. Briefly, pathological findings from neck dissections are treated as the gold standard, overriding any conflicting clinical diagnoses. For levels not dissected, PET-CT is typically the primary source for determining the most likely state of involvement.


## Data Availability

The complete dataset, including additional patients with tumors in primary sites other than the oropharynx, is publicly accessible. It can be [downloaded from LyProX](https://lyprox.org/patients/dataset), where it is also available for interactive exploration, [from GitHub](https://github.com/rmnldwg/lydata), or [from Zenodo](https://zenodo.org/search?q=lydata). Additional details about the datasets and data format are provided in designated Data-in-Brief publications [@ludwig_dataset_2022;@ludwig_multicentric_2024]. These publications do not include the latest data from HVH, which will be described in a future publication.


## Patterns of Contralateral Involvement {#sec-data-strat}

The datasets enable analysis of correlations between contralateral LNL involvement and key risk factors. In @fig-data-strat, we illustrate the prevalence of contralateral LNL involvement, stratified by T-category, the number of ipsilaterally involved LNLs, and whether the tumor extends across the mid-sagittal plane.

![Contralateral involvement stratified by T-category (top left panel), the number of metastatic LNLs ipsilaterally (top right panel), and whether the primary tumor extended over the mid-sagittal plane or was clearly lateralized (bottom left panel). In the bottom right panel, we consider lateralized tumors only, and compare the contralateral involvement prevalence for selected scenarios that vary in their T-category and ipsilateral involvement extent.](figures/fig_data_strat.svg){#fig-data-strat}

### Midline Extension

The bottom left panel of @fig-data-strat shows that tumors crossing the mid-sagittal plane have a substantially higher prevalence of contralateral involvement compared to clearly lateralized tumors. This aligns with the anatomy of the head and neck lymphatic system, which is symmetric, with no major lymph vessels crossing the midline. Interstitial fluids from the primary tumor, presumed to carry malignant cells, can merely diffuse to contralateral lymphatic vessels over short distances. Thus, contralateral spread is more likely when the tumor approaches or crosses the midline. <!-- citation(s) needed? or is this intuition good enough? -->

### T-Category

The top left panel shows a correlation between T-category and contralateral involvement, reflecting T-category's role as a surrogate for the time elapsed between disease onset and diagnosis. Advanced T-category tumors (e.g., T4) generally represent a longer disease progression timeline, providing more opportunity for metastatic spread compared to smaller tumors (e.g., T1).

### Ipsilateral Involvement

The top right panel reveals a positive correlation between ipsilateral and contralateral metastases. Extensive ipsilateral involvement likely indicates a longer or faster disease progression. Additionally, it has been hypothesized that bulky ipsilateral nodal disease may reroute lymphatic drainage toward the contralateral side, potentially increasing the probability of contralateral metastasis. <!-- citation needed? -->

### Correlation of Risk Factors

Midline extension, T-category, and ipsilateral involvement are interrelated risk factors for contralateral metastasis. For instance, {{< var data.late_with_midext_percent >}}% of advanced T-category tumors exhibit midline extension compared to {{< var data.early_with_midext_percent >}}% of early T-category tumors. While the higher fraction of midline extensions in advanced T-category patients partially explains the higher contralateral metastasis rates, T-category itself and ipsilateral involvement also play an additional role.

The bottom right panel of @fig-data-strat considers only patients with lateralized tumors that do not cross the midline. Among early T-category patients with no ipsilateral nodal involvement (levels I-V), only {{< var data.early_lateral_ipsin0_cII_percent >}}% ({{< var data.early_lateral_ipsin0_cII_match >}} of {{< var data.early_lateral_ipsin0_cII_total >}} patients) show involvement in contralateral level II. This proportion increases to {{< var data.early_lateral_ipsiII_cII_percent >}}% ({{< var data.early_lateral_ipsiII_cII_match >}} of {{< var data.early_lateral_ipsiII_cII_total >}}) if ipsilateral level II is involved, to {{< var data.early_lateral_ipsiIIandIII_cII_percent >}}% ({{< var data.early_lateral_ipsiIIandIII_cII_match >}} of {{< var data.early_lateral_ipsiIIandIII_cII_total >}}) if ipsilateral levels II and III are involved, and further to {{< var data.late_lateral_ipsiIIandIII_cII_percent >}}% ({{< var data.late_lateral_ipsiIIandIII_cII_match >}} of {{< var data.late_lateral_ipsiIIandIII_cII_total >}}) for advanced T-category tumors with ipsilateral levels II and III involved.

## Requirements for a Bilateral Model {#sec-requirements}

Based on the observations in @sec-data-strat above, any model predicting the risk of contralateral nodal involvement should account for the following:

1. **Midline Extension:** Tumors extending across the mid-sagittal plane should result in a significantly higher probability of contralateral metastases.
2. **T-Category:** Advanced T-category should correspond to an increased risk of nodal disease. In the hidden Markov model this can be modeled using the expected time of diagnosis, as demonstrated previously for the ipsilateral model [@ludwig_hidden_2021].
3. **Ipsilateral Involvement:** The model should be able to capture the correlation between the extent of ipsi- and contralateral involvement. I.e., a more severe ipsilateral involvement should indicate a higher risk for contralateral metastases.

# Unilateral Model for Lymphatic Progression {#sec-unilateral}

This paper builds on the previously developed unilateral model for ipsilateral lymph node involvement presented in [@ludwig_modelling_2024]. Below, we briefly recap the unilateral model to introduce the notation required for extending the framework to a bilateral model, described in @sec-ext-to-contra. For further details on the ipsilateral model, we refer to earlier publications [@ludwig_hidden_2021;@ludwig_modelling_2024].

We represent a patient’s state of involvement at an abstract time-step $t$ as a vector of hidden binary random variables, where each component corresponds to a lymph node level (LNL):

$$
\mathbf{X}[t] = \begin{pmatrix} X_v[t] \end{pmatrix} \qquad v \in \left\{ 1, 2, \ldots, V \right\}
$$ {#eq-state-def}

Here, $V$ is the number of LNLs the model considers. The values a LNL may take on are $X_v[t] = 0$ (`False`), meaning the LNL $v$ is free of metastatic disease, or $X_v[t] = 1$ (`True`), corresponding to the presence of clinically detected metastases (i.e., occult or macroscopic disease). In total, there are $2^V$ distinct possible lymphatic involvement patterns, which we enumerate from $\boldsymbol{\xi}_0 = \begin{pmatrix} 0 & 0 & \cdots & 0 \end{pmatrix}$ to $\boldsymbol{\xi}_{2^V} = \begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}$. Each LNL's state is observed via another binary random variable $Z_v$ that describes the clinical involvement of a LNL based on imaging: $Z_v = 0$ (`False`) indicates that the LNL $v$ is healthy based on clinical diagnosis, and $Z_v = 1$ (`True`) indicates that LNL $v$ was classified as involved. $X_v$ and $Z_v$ are connected through the sensitivity and specificity of the diagnositc modality.

Based on this, our HMM is fully described by defining the following three quantities:

1. A starting state $\mathbf{X}[t=0]$ at time $t=0$ just before the patient's tumor formed. In our case, this is always the state $\boldsymbol{\xi}_0$ where all LNLs are still healthy.
2. The *transition matrix*
   $$
   \mathbf{A} = \left( A_{ij} \right) = \big( P \left( \mathbf{X}[t+1] = \boldsymbol{\xi}_j \mid \mathbf{X}[t] = \boldsymbol{\xi}_i \right) \big)
   $$ {#eq-trans-matrix}
   where the value at row $i$ and column $j$ represents the probability to transition from state $\boldsymbol{\xi}_i$ to $\boldsymbol{\xi}_j$ during the time-step from $t$ to $t+1$. Note that we prohibit self-healing, meaning that during a transition, no LNL may change their state from $X_v[t]=1$ to $X_v[t+1]=0$. Consequently, many elements of the transition matrix are zero.
3. Lastly, the *observation matrix*
   $$
   \mathbf{B} = \left( B_{ij} \right) = \big( P \left( \mathbf{Z} = \boldsymbol{\zeta}_j \mid \mathbf{X}[t_D] = \boldsymbol{\xi}_i \right) \big)
   $$ {#eq-obs-matrix}
   where in row $i$ and at column $j$ we find the probability to *observe* a lymphatic involvement pattern $\mathbf{Z} = \boldsymbol{\zeta}_j$, given that the true (but hidden) state of involvement at the time of diagnosis $t_D$ is $\mathbf{X}[t_D] = \boldsymbol{\xi}_i$.

The transition matrix $\mathbf{A}$ is parameterized using a directed acyclic graph (DAG) that represents the underlying lymphatic network. Edges from the primary tumor to an LNL are associated with a probability $b_v$ for direct spread to LNL $v$ during one time step. Arcs from a LNL $v$ to a LNL $r$ are parameterized with the probability rate $t_{vr}$ representing the probability of spread to an LNL $r$ that receives efferent lymphatic spread from LNL $v$. In this paper, we build on the DAG shown in @fig-full-graph which was obtained by maximizing the model evidence as described in [@ludwig_modelling_2024].

![Directed acyclic graph (DAG) representing the abstract lymphatic network in the head and neck region. Blue nodes are the LNLs' hidden random variables, the red node represents the tumor, and the orange square nodes depict the binary observed variables. Red and blue arcs symbolize the probability of lymphatic spread along that edge during one time-step. The orange arcs represent the sensitivity and specificity of the observational modality (e.g. CT, MRI, pathology, ...).](static/full-graph.svg){#fig-full-graph width=40%}

Let us now consider the probability distribution over all possible hidden states $\mathbf{X}[t]$ at time $t$. We can get to this distribution by evolving the healthy starting state $\mathbf{X}[t=0] = \boldsymbol{\xi}_0$ at time $t=0$ step by step, by successively multiplying this vector with the transition matrix $\mathbf{A}$: $\mathbf{X}[t+1] = \mathbf{X}[t] \cdot \mathbf{A}$. For later use, we define at this point a matrix $\boldsymbol{\Lambda}$ that collects these distributions for all considered times:

$$
\boldsymbol{\Lambda} = P \left( \mathbf{X} \mid \mathbf{t} \right) = \begin{pmatrix}
\boldsymbol{\pi}^\intercal \cdot \mathbf{A}^0 \\
\boldsymbol{\pi}^\intercal \cdot \mathbf{A}^1 \\
\vdots \\
\boldsymbol{\pi}^\intercal \cdot \mathbf{A}^{t_\text{max}} \\
\end{pmatrix}
$$ {#eq-lambda-matrix}

where the $k$-th row in this matrix corresponds to the probability distribution over hidden states after $t=k-1$ time-steps.

At the time of diagnosis, $0 \leq t_D \leq t_\text{max}$, we multiply the evolved state distribution with the observation matrix $\mathbf{B}$ to obtain the distribution over all possible diagnoses. However, the exact time of diagnosis, $t_D$, is unknown; that is, we do not know the number of time-steps over which the HMM should be evolved. To address this, we marginalize over all possible diagnosis times, allowing the diagnosis to occur at any time-step, albeit with different weights. These weights are defined by a prior distribution over $t_D$, which can vary depending on the patient's T-category. For example, the time-prior for early T-category patients, $P(t_D \mid \text{early})$, may put more weight on earlier time-steps, reflecting -- on average -- earlier detection, compared to the prior for advanced T-category patients, $P(t_D \mid \text{advanced})$.

The probability distribution over $\mathbf{X}$ for the example of an early T-category patient is given by

$$
P\left( \mathbf{X} \mid \text{T}x = \text{early} \right) = \sum_{t=0}^{t_\text{max}} P \left( \mathbf{X} \mid t \right) \cdot P(t \mid \text{T}x = \text{early})
$$

In this work, we use binomial distributions $\mathfrak{B} \left( t_D, p_{\text{T}x} \right)$ as time-priors which have one free parameter $p_{\text{T}x}$ for each group of patients we differentiate based on T-category. Also, we fix $t_\text{max} = 10$, which means that the expected number of time-steps from the onset of a patient's disease to their diagnosis is $\mathbb{E}\left[ t_D \right] = 10 \cdot p_{\text{T}x}$.


## Likelihood Function of the Unilateral Model

The probabiltiy for a patient to present with a diagnosis $\mathbf{Z} = \boldsymbol{\zeta}_i$ a T-category $\text{T}x$ tumor can now be written as:

$$
\ell = P \left( \mathbf{Z} = \boldsymbol{\zeta}_i \mid \text{T}x \right) = \sum_{t=0}^{t_\text{max}} \left[ \boldsymbol{\xi}_0 \cdot \mathbf{A}^t \cdot \mathbf{B} \right]_i \cdot P \left( t \mid \text{T}x \right)
$$ {#eq-single-patient-llh}

here, $\left[ \ldots \right]_i$ we denote the $i$-th component of the vector in the square brackets. Note that it is also possible to account for missing involvement information: If a diagnosis (like fine needle aspiration (FNA)) is only available for a subset of all LNLs, we can sum over all those possible complete observed states $\boldsymbol{\zeta}_j$ that match the provided diagnosis.

The term above represents a single patient's contribution to the overall likelihood function that is a product of such terms for each patient. This single-patient likelihood $\ell$ in @eq-single-patient-llh depends on the spread parameters shown in @fig-full-graph via the transition matrix $\mathbf{A}$ and on the binomial parameters $p_{\text{T}x}$ via time-priors. In this work, we will only differentiate between "early" (T1 & T2) and "advanced" (T3 & T4) T-categories. Therefore, the parameter space of the unilateral model is:

$$
\boldsymbol{\theta} = \left( \left\{ b_v \right\}, \left\{ t_{vr} \right\}, p_\text{early}, p_\text{adv.} \right) \quad \text{with} \quad \genfrac{}{}{0pt}{2}{v\leq V}{r\in\operatorname{pa}(v)}
$$ {#eq-param-space}

And it is our goal to infer optimal parameter values of from a given dataset $\mathcal{D}$ (consisting of diagnoses and T-categories) of OPSCC patients. The likelihood to observe this cohort of $N$ patients, given a set of parameters $\boldsymbol{\theta}$ is simply the product of their individual likelihoods as defined in @eq-single-patient-llh. For numerical reasons, we typically compute the data likelihood in log space:

$$
\log \mathcal{L} \left( \mathcal{D} \mid \boldsymbol{\theta} \right) = \sum_{i=1}^N \log \ell_i
$$ {#eq-log-likelihood}

The methodology we use to infer the model's parameters is detailed in @sec-sampling.


# Extension to a Bilateral Model {#sec-ext-to-contra}

A straightforward approach to modeling contralateral lymphatic spread would be to use two independent unilateral models, as described in @sec-unilateral, possibly with shared parameters such as the distribution of diagnosis times or spread between LNLs ($t_{vr}$). However, this method would fail to capture the correlation between ipsilateral and contralateral involvement discussed in @sec-data-strat, particularly the observed increase in contralateral involvement with greater severity of ipsilateral spread.

Thus, we extend the formalism in @sec-unilateral in such a way that the model's ipsi- and contralateral side evolve synchronously over time. To achieve that, we start by writing down the posterior distribution of involvement, which is now a joint probability of an involvement $\mathbf{X}^\text{i}$ ipsilaterally *and* an involvement $\mathbf{X}^\text{c}$ contralaterally, given a diagnosis of the ipsilateral LNLs $\mathbf{Z}^\text{i}$ and of the contralateral ones $\mathbf{Z}^\text{c}$:

$$
P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \mid \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \right) = \frac{P \left( \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \mid \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right) P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right)}{P \left( \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \right)}
$$ {#eq-bilateral-bayes}

For the sake of brevity, we omit the dependency on the parameters and the T-category here.

The probability of the diagnoses given a hidden state factorises: $P \left( \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \mid \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right) = P \left( \mathbf{Z}^\text{i} \mid \mathbf{X}^\text{i} \right) \cdot P \left( \mathbf{Z}^\text{c} \mid \mathbf{X}^\text{c} \right)$, and the two factors are described through observation matrices $\mathbf{B}^\text{i}$ and $\mathbf{B}^\text{c}$.

The term representing the model's prior probability of hidden involvement does not factorize. However, we assume no direct lymphatic drainage from ipsilateral to contralateral LNLs, as major lymph vessels do not cross the mid-sagittal plane. In the graphical model, this translates to the absence of directed arcs between ipsilateral and contralateral LNLs, implying that contralateral tumor spread occurs solely via the primary tumor. We can thus write the joint probability $P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right)$ as a factorising sum:

$$
\begin{aligned}
P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right) &= \sum_{t=0}^{t_\text{max}} P(t) \cdot P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \mid t \right) \\
&= \sum_{t=0}^{t_\text{max}} P(t) \cdot P \left( \mathbf{X}^\text{i} \mid t \right) \cdot P \left( \mathbf{X}^\text{c} \mid t \right)
\end{aligned}
$$ {#eq-bilateral-marginal}

This assumption is intuitive: since no major lymph vessels cross the midline, the ipsilateral and contralateral sides of the lymphatic network evolve independently over time. However, they are indirectly coupled through time. For example, a joint state with severe contralateral involvement and limited ipsilateral involvement is improbable: Severe contralateral involvement typically occurs at later time steps, when limited ipsilateral involvement is unlikely.

Using @eq-bilateral-marginal along with @eq-lambda-matrix, we can write the above distribution algebraically as a product:

$$
P \left( \mathbf{X}^\text{i} = \boldsymbol{\xi}_n, \mathbf{X}^\text{c} = \boldsymbol{\xi}_m \right) = \left[ \boldsymbol{\Lambda}^\intercal_\text{i} \cdot \operatorname{diag} P(\mathbf{t}) \cdot \boldsymbol{\Lambda}_\text{c} \right]_{n,m}
$$ {#eq-bilateral-marginal-algebra}


## Parameter Symmetries {#sec-params-symmetry}

The matrices $\boldsymbol{\Lambda}_\text{i}$ and $\boldsymbol{\Lambda}_\text{c}$ could, in principle, be parameterized with entirely separate parameters, allowing ipsilateral and contralateral spread rates to differ substantially. However, we simplify the parameter space by sharing parameters between the two sides, based on the following three assumptions:

1. **Shared Graph Structure**: Both ipsilateral and contralateral spread are described by the same graph shown in @fig-full-graph.
2. **Symmetric Spread Among LNLs**: The spread among LNLs is assumed to be the same on both sides, reflecting the symmetric structure of the lymphatic system. Consequently, the spread rates between nodes should also be symmetric. This is formalized as:\
  $$
  t_{rv}^\text{c} = t_{rv}^\text{i}
  $$ {#eq-symmetries}
  for all $v \leq V$ and $r \in \operatorname{pa}(v)$ being all nodes that spread to $v$.
1. **Asymmetric Spread from Tumor**: The probabilities of direct spread from the primary tumor are clearly different for the ipsi- and contralateral neck. In addition, tumor spread to the contralateral side varies depending on whether the tumor crosses the mid-sagittal plane. This would result in three sets of rates for tumor spread to the LNLs: (1) the spread to ipsilateral LNLs $b^\text{i}_v$, (2) the spread to contralateral LNLs as long as the tumor is lateralized $b_v^{\text{c},\epsilon=\texttt{False}}$, (3) the spread to contralateral LNLs when the tumor crosses the midline $b_v^{\text{c},\epsilon=\texttt{True}}$. In this work, however, we chose to define the latter set as a linear mix of ipsilateral tumor spread and contralateral spread in case of a clearly lateralized tumor. Thus, using $\alpha \in [0,1]$ as this mixing parameter, we have
  $$
  b_v^{\text{c},\epsilon=\texttt{True}} = \alpha \cdot b_v^\text{i} + (1 - \alpha) \cdot b_v^{\text{c},\epsilon=\texttt{False}}
  $$ {#eq-mixing}
  In this way, the tumor's midline extension causes the contralateral spread to become more like the spread to the ipsilateral side.

The full parameter space of this model is now:

$$
\boldsymbol{\theta} = \left( \left\{ b_v^\text{i} \right\}, \left\{ b_v^\text{c} \right\}, \alpha, \left\{ t_{vr} \right\}, p_\text{early}, p_\text{adv.} \right) \quad \text{with} \quad \genfrac{}{}{0pt}{2}{v\leq V}{r\in\operatorname{pa}(v)}
$$ {#eq-bi-param-space}

This results in less than a doubling of parameters compared to the unilateral model. From these parameters, we construct three transition matrices: the unchanged $\mathbf{A}_\text{i}$ for the ipsilateral side, $\mathbf{A}_\text{c}^{\epsilon=\texttt{False}}$ for contralateral progression while the tumor is lateralized, and $\mathbf{A}_\text{c}^{\epsilon=\texttt{True}}$ for cases where the tumor crosses the mid-sagittal plane.


## Modelling Midline Extension {#sec-midline}

Most tumors crossing the midline at the time of diagnosis likely began as lateralized tumors that grew over the midline at a later point in time. As a result, the transition matrix $\mathbf{A}_\text{c}^{\epsilon=\texttt{True}}$ applies only to a subset of time-steps.

To account for this, we model the tumor's extension over the mid-sagittal plane as an additional binary random variable $\epsilon$. A tumor starts as lateralized, with a finite probability $p_\epsilon$ at each time step of crossing the midline. The overall probabilities of a patient having a clearly lateralized tumor or one extending over the mid-sagittal plane after $t$ time steps are then given by

$$
\begin{aligned}
P(\epsilon = \texttt{False} \mid t) &= (1 - p_\epsilon)^t \\
P(\epsilon = \texttt{True} \mid t) &= 1 - P(\epsilon = \texttt{False} \mid t)
\end{aligned}
$$

Using this, it is straightforward to write down the matrix of state distributions for all time-steps, as in @eq-lambda-matrix, covering the joint distribution over the contralateral hidden state and the midline extension:

$$
\boldsymbol{\Lambda}_\text{c}^{\epsilon=\texttt{False}} =
\left(
\begin{array}{r}
\boldsymbol{\pi}^\intercal \cdot \left( \mathbf{A}_\text{c}^{\epsilon=\texttt{False}} \right)^{0\phantom{t_\text{max}}} \\
(1-p_\epsilon) \cdot \boldsymbol{\pi}^\intercal \cdot \left( \mathbf{A}_\text{c}^{\epsilon=\texttt{False}} \right)^{1\phantom{t_\text{max}}} \\
\hfill \vdots \hfill \\
(1-p_\epsilon)^{t_\text{max}} \cdot \boldsymbol{\pi}^\intercal \cdot \left( \mathbf{A}_\text{c}^{\epsilon=\texttt{False}} \right)^{t_\text{max}\phantom{0}} \\
\end{array}
\right)
$$

here, we used the transition matrix $\mathbf{A}_\text{c}^{\epsilon=\texttt{False}}$ that depends on the base spread parameters $b_v^{\text{c},\epsilon=\texttt{False}}$.

The case of midline extension is more complex: we already marginalize over the exact time step when the tumor grows over the mid-sagittal plane. However, at the point of crossing, the contralateral transition matrix must switch to the increased spread rates, $b_v^{\text{c}, \epsilon=\texttt{True}}$, as defined by the linear mixing in @eq-mixing. To correctly perform this marginalization, we iteratively construct the joint distribution $P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid t \right)$.

We begin at $t=0$, where all contralateral LNLs are healthy (i.e., $\mathbf{X}_\text{c}=\boldsymbol{\xi}_0$) and the tumor is lateralized ($\epsilon=\texttt{False}$):

$$
P \left( \mathbf{X}^\text{c} = \boldsymbol{\xi}_0, \epsilon=\texttt{False} \mid t=0 \right) = 1
$$

while all other states have zero probability.

At some later time step $t=\tau+1$, there are two scenarios to marginalize over:

1. **The tumor was lateralized at $t=\tau$ and grew over the midline at $t=\tau+1$:**  
   In this case, the probability of midline extension at $t=\tau+1$ is $p_\epsilon$. This probability weights the contralateral state distribution that had previously evolved without increased contralateral spread.

2. **The tumor had already crossed the midline before $t=\tau$:**  
   Here, the tumor remains in the midline-crossed state with probability 1. To account for this scenario, we simply include the distribution $P\left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid \tau \right)$ from the previous time step.

Combining these scenarios leads to a recursive formulation:

$$
\begin{multlined}
P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid \tau + 1 \right) \\
= \big[ p_\epsilon P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{False} \mid \tau \right) + P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid \tau \right) \big]^\top \cdot \mathbf{A}_\text{c}^{\epsilon=\texttt{True}}
\end{multlined}
$$

We can collect the iteratively computed distributions for the midline extension case to define the matrix over the states given all time-steps, as in @eq-lambda-matrix:

$$
\boldsymbol{\Lambda}_\text{c}^{\epsilon=\texttt{True}} = \begin{pmatrix}
P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid 0 \right) \\
P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid 1 \right) \\
\vdots \\
P \left( \mathbf{X}^\text{c}, \epsilon=\texttt{True} \mid t_\text{max} \right) \\
\end{pmatrix}
$$

In analogy to @eq-bilateral-marginal-algebra, we can now write the joint distribution of ipsi- and contralateral involvement and midline extension algebraically:

$$
P \left( \mathbf{X}^\text{i} = \boldsymbol{\xi}_n, \mathbf{X}^\text{c} = \boldsymbol{\xi}_m, \epsilon \right) = \left[ \boldsymbol{\Lambda}^\intercal_\text{i} \cdot \operatorname{diag} P(\mathbf{t}) \cdot \boldsymbol{\Lambda}_\text{c}^\epsilon \right]_{n,m}
$$ {#eq-midline-marginal-algebra}

With the above, we compute the likelihood of all patients with and without midline extension separately. And if for some patients the information of tumor lateralization is not available, we can simply marginalize over the unknown variable $\epsilon \in \{ \texttt{False}, \texttt{True} \}$.

The final parameter space of our extended model has now reached this size:

$$
\boldsymbol{\theta} = \left( \left\{ b_v^\text{i} \right\}, \left\{ b_v^\text{c} \right\}, \alpha, \left\{ t_{vr} \right\}, p_\text{early}, p_\text{adv.}, p_\epsilon \right) \quad \text{with} \quad \genfrac{}{}{0pt}{2}{v\leq V}{r\in\operatorname{pa}(v)}
$$ {#eq-ext-param-space}


## Model Prediction in the Bayesian Context

Our stated goal is to compute the risk for a patient's true ipsi- and contralateral nodal involvement states $\mathbf{X}^\text{i}$ and $\mathbf{X}^\text{c}$, *given* their individual diagnosis $d = \left( \boldsymbol{\zeta}^\text{i}_k, \boldsymbol{\zeta}^\text{c}_\ell, \epsilon, \text{T}x \right)$. Here, this diagnosis consists of the observed ipsi- and contralateral nodal involvements, the patient's midline extension $\epsilon$, and their tumor's T-category $\text{T}x$. Using Bayes' law, we can write this risk as:

$$
P \big( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \mid d, \boldsymbol{\hat{\theta}} \big)
= \frac{P \left( \boldsymbol{\zeta}^\text{i}_k \mid \mathbf{X}^\text{i} \right) P \left( \boldsymbol{\zeta}^\text{c}_\ell \mid \mathbf{X}^\text{c} \right) P \big( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \mid \boldsymbol{\hat{\theta}}, \text{T}x \big)}
{\sum_{i=0}^{2^V} \sum_{j=0}^{2^V} \mathcal{C}_{ij}}
$$ {#eq-uni-bayes-law}

with the normalization constants

$$
\mathcal{C}_{ij} = P \left( \boldsymbol{\zeta}^\text{i}_k \mid \mathbf{X}^\text{i}=\boldsymbol{\xi}^\text{i}_i \right) P \big( \boldsymbol{\zeta}^\text{c}_\ell \mid \mathbf{X}^\text{c}=\boldsymbol{\xi}^\text{c}_j \big) P \big( \mathbf{X}^\text{i}=\boldsymbol{\xi}^\text{i}_i, \mathbf{X}^\text{c}=\boldsymbol{\xi}^\text{c}_j, \epsilon \mid \boldsymbol{\hat{\theta}}, \text{T}x \big)
$$

The terms $P \left( \boldsymbol{\zeta}^\text{i}_k \mid \mathbf{X}^\text{i} \right)$ and $P \left( \boldsymbol{\zeta}^\text{c}_\ell \mid \mathbf{X}^\text{c} \right)$ are defined solely by sensitivity and specificity of the diagnostic modality. These terms already appeared in the definition of the observation matrx in @eq-obs-matrix. The *prior* $P \big( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \mid \boldsymbol{\hat{\theta}}, \text{T}x \big)$ in the above equation is the crucial term that is supplied by a trained model and its parameters $\boldsymbol{\hat{\theta}}$.

It is possible to compute this *posterior* probability of true involvement not only for one fully defined state $(\mathbf{X}^\text{i}, \mathbf{X}^\text{c})$, but also for e.g. individual LNLs: For example, the risk for involvement in the contralateral level IV would be a marginalization over all combination of ipsi- states $\boldsymbol{\xi}^\text{i}_i$ contralateral states $\boldsymbol{\xi}^\text{c}_j$ where $\xi^\text{c}_{j4}=1$. Formally:

$$
\begin{multlined}
P \big( \text{IV}^\text{c} \mid \mathbf{Z}^\text{i}=\boldsymbol{\zeta}^\text{i}_k, \mathbf{Z}^\text{c}=\boldsymbol{\zeta}^\text{c}_\ell, \boldsymbol{\hat{\theta}}, \text{T}x \big) \\
= \sum_k \sum_{\ell \, : \, \xi_{\ell 4}=1} P \big( \mathbf{X}^\text{i} = \boldsymbol{\xi}^\text{i}_k, \mathbf{X}^\text{c} = \boldsymbol{\xi}^\text{c}_\ell \mid \boldsymbol{\zeta}^\text{i}_k, \boldsymbol{\zeta}^\text{c}_\ell, \epsilon, \boldsymbol{\hat{\theta}}, \text{T}x \big)
\end{multlined}
$$ {#eq-marg-over-posterior}


# Computational Methods {#sec-methods}

This section details the experimental setup. All figures, tables, and results are fully reproducible via the GitHub repository [rmnldwg/bilateral-paper].

[rmnldwg/bilateral-paper]: https://github.com/rmnldwg/bilateral-paper

## Training Data {#sec-training-data}

We trained the model using the dataset of {{< var data.num_patients >}} patients described in @sec-data. The consensus decision on lymphatic involvement is assumed to correspond to the true hidden state of involvement $\mathbf{X}$. Patients with T1 and T2 category tumors have been grouped into an "early" T-category group, those with T3 and T4 tumors into the "advanced" T-category group.

## MCMC Sampling {#sec-sampling}

We used the Python package [`emcee`] [@foreman-mackey_emcee_2013] for parameter inference, implementing efficient MCMC sampling with parallel affine-invariant samplers. The sampling algorithms employed differential evolution moves [@terbraak_differential_2008;@nelson_run_2013], with the likelihood implemented by our [`lymph-model`] Python package.

We initialized {{< var sampling.walkers_per_dim >}} parallel samplers ("walkers") with random values from the unit cube, effectively representing a uniform prior distribution over the model parameters. Convergence was determined by two criteria:

1. The change in autocorrelation time was less than {{< var sampling.rel_thresh >}}.
2. The autocorrelation estimate dropped below $n$ / {{< var sampling.trust_fac >}}, where $n$ is the chain length. Earlier autocorrelation estimates might not be trustworthy.

Samples from this *burn-in phase* before convergence were discarded. After that, we drew {{< var sampling.nsteps >}} additional samples, spaced {{< var sampling.thin >}} steps apart.

We verified sampling convergence in @fig-model-burnin-history by examining the MCMC chain's autocorrelation time and walker acceptance fractions.

[`emcee`]: https://emcee.readthedocs.io/en/stable/
[`lymph-model`]: https://lymph-model.readthedocs.io/en/stable/

![Burn-in phase monitoring of MCMC sampling. Left: Estimated autocorrelation time, indicating converging when stable and below the trust threshold. Right: Average acceptance fraction of parallel walkers, with ~30% indicating good mixing.](figures/fig_model_burnin_history.svg){#fig-model-burnin-history}


## Comparing the Observed and Predicted Prevalence of Involvement Patterns {#sec-prevalence}

We evaluate the model's ability to describe the observed frequencies of lymphatic involvement patterns. We compare the prevalence of selected involvement patterns in the data to the model's predicted prevalence, given patient scenarios. A "scenario" includes the patient's T-category $\text{T}x$ and whether the tumor extended over the mid-sagittal plane, i.e. $\epsilon=\texttt{True}$ or $\epsilon=\texttt{False}$. An involvement pattern specifies all ipsi- and contralateral LNLs' status as "healthy", "involved", or "masked" (ignored). 

For example, in @fig-model-prevalences-overall (top left panel) we assess contralateral LNL II involvement prevalence for early T-category (T0-T2) and no midline extension ($\epsilon=\texttt{False}$). In the data, $n=$ {{< var data.early_lateral_cII_total >}} such patients were observed, with $k=$ {{< var data.early_lateral_cII_match >}} exhibiting contralateral LNL II involvement -- an observed prevalence of $q=$ {{< var data.early_lateral_cII_percent >}}%. To visualize the data prevalence, we plot a *beta posterior* over $q$ -- with a uniform beta prior -- multiplied with the binomial likelihood for $k$ out of $n$ patients, given $q$. The resulting distribution has its maximum at $q=k / n$ and nicely captures the statistical uncertainty in the observed cohort: In the top right panel of @fig-model-prevalences-overall, we consider the case of early T-category tumors extending over the midline. The dataset contains only {{< var data.early_with_midext_cII_total >}} such patients, out of which {{< var data.early_with_midext_cII_match >}} had contralateral LNL II involvement. Consequently, the beta distribution over the observed prevalence is much wider.

The model's predicted prevalence to compare it with is computed as:

$$
\begin{multlined}
P \left( \text{II}^\text{c} \mid \epsilon=\texttt{False}, \text{T}x=\text{early} \right) \\ = \frac{ \sum_k \sum_{\ell \, : \, \xi_{\ell 2}=1}  P \left( \mathbf{X}^\text{i}=\boldsymbol{\xi}_k^\text{i}, \mathbf{X}^\text{c}=\boldsymbol{\xi}_\ell^\text{c}, \epsilon=\texttt{False} \mid \text{T}x=\text{early} \right)}{ \sum_k \sum_\ell P \left(\mathbf{X}^\text{i}=\boldsymbol{\xi}_k^\text{i}, \mathbf{X}^\text{c}=\boldsymbol{\xi}_\ell^\text{c}, \epsilon=\texttt{False} \mid \text{T}x=\text{early} \right)}
\end{multlined}
$$

In the enumator, we marginalize over all combinations of states in both sides of the neck where the contralateral LNL II is involved. This is similar to the marginalization in @eq-marg-over-posterior, although we are summing over different quantities. In the denominator, we simply sum out all LNL involvement, leaving only the joint distribution over midline extension and diagnose time $P \left( \epsilon, t \right)$ marginalized over $t$ using the early T-category's time-prior.

We display model predictions as histograms, each value computed from one of the MCMC samples. Ideally, these approximate the location and width of the beta posteriors from the data showing an accurate and precise fit.

Note that we omit the y-axis in these figures, as their numerical value is not intuitively interpretable. We instead use the free space to label e.g. rows in an array of subplots.


# Results: Model evaluation {#sec-results}

In @tbl-midline-params, we tabulate the mean and standard deviation of the sampled parameters. The bilateral model mostly reproduces the ipsilateral spread parameter values reported in the earlier publication on the unilateral model [@ludwig_modelling_2023]. Any discrepancies may be due to differences in the patient cohorts. Therefore, we omit the analysis of the ipsilateral spread patterns and focus instead on the analysis of contralateral involvement. The small contralateral spread parameters $b^c_v$ compared to the ipsilateral parameters $b^i_v$ adequately reflect the low prevalence of contralateral lymph node involvement for lateralized tumors. The mixing parameter $\alpha$ ({{< var model.params_percent.mixing >}}%) describes that the probability of contralateral spread is higher for tumors extending over the mid-sagittal plane, but still lower than the probability for ipsilateral spread.

The high value of $t_{23}$ ({{< var model.params_percent.IItoIII_spread >}}%) compared to the value of $b_3^\text{i}$ ({{< var model.params_percent.ipsi_TtoIII_spread >}}%) shows that the ipsilateral LNL III is rarely involved without the upstream involvement of LNL II.

```{python}
#| echo: false
#| label: tbl-midline-params
#| tbl-cap: Mean sampled parameter estimates of the midline model and the respective standard deviation. The parameters set to fixed values are the maximum number of time steps $t_{max}=10$ and the time prior parameter for early T-category patients $p_{early}=0.3$.
def map_param_names(name: str) -> str:
  """Make parameter names more readable."""
  try:
    return {
      "midext_prob": "Mid. ext. probability",
      "mixing": "Mixing ⍺",
      "late_p": "late T-cat. binom. prob.",
    }[name]
  except KeyError:
    return re.sub(
      r"(ipsi|contra)_", r"\g<1>: ", name
    ).replace(
      "to", " ➜ "
    ).replace(
      "_spread", ""
    )

model = shared.get_model(which="full", load_samples=True)
samples = shared.get_samples(which="full")

names = [map_param_names(p) for p in model.get_params()]
means, stds = samples.mean(axis=0), samples.std(axis=0)

early_midext_prob = model.state_dist(t_stage="early")[1].sum()
late_midext_prob = model.state_dist(t_stage="late")[1].sum()

params_table = pd.DataFrame({"Parameter": names, "Mean": means, "Std. Dev.": stds})
(
  params_table.style
  .format("{:.2%}", subset=["Mean"])
  .format("± {:.2%}", subset=["Std. Dev."])
  .apply(right_align, subset=["Mean", "Std. Dev."])
  .hide()
)
```

## Illustration of the model

In this subsection, we illustrate key aspects of the mathematical framework introduced earlier. The top panel of @fig-model-midext-evo shows the prior distribution over diagnosis times, $P(t)$. Based on the parameterization, early T-category tumors are on average diagnosed after {{< var model.early_expected_time >}} time steps, while advanced T-category tumors are diagnosed later, averaging {{< var model.late_expected_time >}} time steps. This is due to the learned value of $p_\text{adv.}$, which is {{< var model.params_percent.late_p >}}% for advanced T-category tumors. The tumor's average probability per time step of growing over the midline, $p_\epsilon$, was found to be {{< var model.params_percent.midext_prob >}}%. Using this value, the conditional probability of midline extension, $P(\epsilon \mid t)$, can be computed for a given time step $t$ (red line in the top panel of @fig-model-midext-evo). The bottom panel visualizes the joint probability $P(\epsilon, t)$, showing the likelihood of diagnosis at time $t$ with specific states of midline extension and T-category.

![The top panel shows the prior probability to be diagnosed at time step $t$ for early and late T-category tumors as bars. The conditional probability of midline extension ($\epsilon=\texttt{True}$) given time step $t$ is shown as a line plot. The bottom panel illustrates the joint probability of being diagnosed at time $t$ and having a tumor that crosses the midline.](figures/fig_model_midext_evo.svg){#fig-model-midext-evo width="50%"}

The framework models the joint probability distribution of midline extension and ipsi- and contralateral lymph node involvement, $P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \right)$. This is visualized in @fig-model-algebra-equation, which represents the calculation defined in @eq-midline-marginal-algebra. To simplify the interpretation, the example focuses only on LNLs II, III, and IV, reducing the state space to $2^3 = 8$ possible states per side, and $2 \times 8 \times 8 = 128$ states in total. LNLs I, V, and VII are excluded, along with their spread parameters, while remaining parameters are set to their mean values from @tbl-midline-params.

![Visual representation of @eq-midline-marginal-algebra. The left and right matrices represent the time evolution of hidden states for the ipsi- and contralateral necks, respectively. The right matrices distinguish between the cases of no midline extension (top) and midline extension (bottom). The central diagonal matrix shows the time-prior for late T-category tumors. This computation yields the joint distribution $P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \right)$, visualized in @fig-model-state-dist.](figures/fig_model_algebra_equation.svg){#fig-model-algebra-equation}

The left matrix in @fig-model-algebra-equation shows the time evolution of the probability distribution over the ipsilateral involvement states, starting from the healthy state $[0,0,0]$. The two right matrices show the contralateral state evolution, distinguishing between the cases of midline extension and no midline extension. At $t=0$, the contralateral neck begins in the healthy state $[0,0,0]$ without midline extension. The central matrix shows the time prior for late T-category tumors. The matrix multiplication results in the joint distribution $P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right)$, visualized in @fig-model-state-dist.

This joint distribution is presented as two heatmaps, corresponding to the two states of midline extension. The most likely state is a lateralized tumor with ipsilateral level II involvement and no contralateral involvement, having a probability of approximately 25%. The next most probable state is a lateralized tumor with ipsilateral levels II and III involved, but without contralateral involvement. The most likely state with contralateral involvement corresponds to tumors with midline extension, showing involvement of contralateral level II and ipsilateral levels II and III.

![The joint distribution $P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \right)$ over ipsi- and contralateral states and midline extension for late T-category tumors. The distribution is shown as two separate heatmaps for the two binary values of the midline extension variable $\epsilon$. These matrices are the result of @eq-midline-marginal-algebra, visualized in @fig-model-algebra-equation.](figures/fig_model_state_dist.svg){#fig-model-state-dist}

## Prevalence predictions for contralateral involvement

The bilateral model was designed to meet the requirements outlined in @sec-requirements. Here, we evaluate the model's ability to quantitatively capture the observed patterns of lymph node involvement in the dataset. Specifically, we compare the model's predictions for contralateral involvement to the observed data across scenarios that vary by T-category, midline extension, and ipsilateral involvement.

### Dependence of Contralateral Involvement on T-Category and Midline Extension

In @fig-model-prevalences-overall, we compare the prevalence of contralateral involvement for LNLs II, III, and IV.

![Comparison of predicted (histograms) vs observed (beta posteriors) prevalences, shown for the contralateral LNLs II (blue), III (orange), and IV (green). The top row shows scenarios with early T-category tumors, the bottom row for late T-category tumors. The left column depicts scenarios where the primary tumor is clearly lateralized, the right column scenarios of tumors extending over the mid-sagittal plane. This figure illustrates the model's ability to describe the prevalence of involvement for different combinations of the risk factors T-category and midline extension.](figures/fig_model_prevalences_overall.svg){#fig-model-prevalences-overall}

@Fig-model-prevalences-overall demonstrates the model's ability to account for the risk factors T-category and midline extension. Consistent with the data, the model predicts that the prevalence of contralateral LNL II involvement increases from {{< var prevalence.early_noext_iN0_cII.predicted >}}% for early T-category lateralized tumors to {{< var prevalence.late_ext_iN0_cII.predicted >}}% for advanced T-category tumors that cross the midline. Similarly, contralateral LNL III involvement rises from around {{< var prevalence.early_noext_iN0_cIII.predicted >}}% for early T-category lateralized tumors to nearly {{< var prevalence.late_ext_iN0_cIII.predicted >}}% for advanced T-category midline-extending tumors.

### Influence of Upstream Involvement on Contralateral Metastasis

![The influence of the upstream LNL II's involvement on the prevalence of contralateral level III for the four combinations of tumor lateralization (lateralized or extending over midline) and T-category (early or advanced). Our model predictions (histograms) are plotted against the observations in the data (beta posteriors).](figures/fig_model_prevalences_upstream.svg){#fig-model-prevalences-upstream}

@Fig-model-prevalences-upstream highlights the influence of upstream LNL II involvement on contralateral LNL III metastasis. The contralateral LNL III rarely harbors metastases when its upstream LNL II is healthy, a correlation well-captured by the model. Out of {{< var data.late_with_midext_upstream_cIII_total >}} patients with an advanced tumor crossing the midline and involvement in the upstream LNL II, {{< var data.late_with_midext_upstream_cIII_match >}} had involvement in the contralateral LNL III. In contrast, out of {{< var data.early_lateral_cIII_total >}} patients with a clearly lateralized early T-category tumor and no upstream involvement, only {{< var data.early_lateral_cIII_match >}} showed contralateral LNL III involvement. This is well captured by our model.

### Correlation between Ipsi- and Contralateral Involvement

![Comparison of the computed and observed prevalences for scenarios that illustrate the model's capability of accounting for the correlation between ipsi- and contralateral involvement. We show three scenarios where we consider the joint involvement of contralateral LNL II together with different ipsilateral involvements: 1) the ipsilateral neck shows no involvement in green (LNLs I to V are healthy, LNL VII is unspecified because data on it is missing for some patients), 2) where ipsilateral LNL II is involved in orange (LNLs I, III, IV, and V are healthy), and 3) where ipsilateral LNLs II and III are involved in red (LNLs I, IV, and V are healthy). These three scenarios are plotted for all combinations of T-category (early in top row, advanced in bottom row) and tumor lateralization (lateralized in left column, extending over mid-sagittal plane in the right column).](figures/fig_model_prevalences_with_ipsi.svg){#fig-model-prevalences-with-ipsi}

In @fig-model-prevalences-with-ipsi, the model's ability to capture the correlation between ipsi- and contralateral involvement is demonstrated. The marginals of the joint distribution highlight contralateral LNL II involvement alongside varying ipsilateral LNL involvement states. Despite having no direct connections between the two sides, the model successfully predicts these correlations, which arise purely through the shared diagnosis time. 

For example, the model accurately predicts that contralateral LNL II involvement is rare when the ipsilateral neck is completely healthy (green histograms). However, if ipsilateral LNL II is involved, contralateral involvement becomes more likely. Notably, the model achieves this without specific parameters to quantify ipsi-contralateral correlations, relying instead on the inherent structure of the time-dependent dynamics.


# Results: Prediction of Risk for Occult Disease {#sec-results-risk}

For clinical applications, the model needs to estimate the risk of occult metastases in clinically negative LNLs based on a patient's diagnosis. The diagnosis includes the T-category, tumor lateralization, and the clinically detected involvement of LNLs based on imaging and possibly fine needle aspiration (FNA). 

We assume imaging detects lymph node involvement with a sensitivity of 81% and a specificity of 76% [@debondt_detection_2007], while FNA has a sensitivity of 80% and a specificity of 98%. This implies that FNA-confirmed involvement is highly reliable, with almost no false positives. 

![Histograms over the predicted risk of occult involvement in contralateral LNL II (top), III (middle), and IV (bottom), shown for various combinations of T-category, tumor lateralization, and clinical LNL diagnoses. All LNLs not explicitly mentioned in the legend, including the LNL for which the risk of occult disease was computed, were assumed to be clinically negative (specificity 76%, sensitivity 81%).](figures/fig_model_risks.svg){#fig-model-risks width="50%"}

## Contralateral LNL II

The predicted risk of occult disease in contralateral LNL II is shown in the left panel of @fig-model-risks. Tumor lateralization is the strongest determinant of risk:  

- A patient with a lateralized early T-category tumor and ipsilateral LNL II involvement has a predicted risk of {{< var risk.early_noext_iII_cN0_II >}}% for occult contralateral LNL II disease (green histogram).  
- For an early T-category tumor that extends over the midline, with ipsilateral LNL II involved, the risk increases to {{< var risk.early_ext_iN0_cN0_II >}}% (orange histogram).  

Advanced T-category further increases risk but has less impact than midline extension. For instance:  

- An early T-category tumor crossing the midline with ipsilateral involvement of LNLs II and III has a {{< var risk.early_ext_iII+III_cN0_II >}}% risk of contralateral LNL II disease (red histogram).  
- For the same scenario but an advanced T-category tumor, the risk rises to {{< var risk.late_ext_iII+III_cN0_II >}}% (purple histogram).  

The degree of ipsilateral involvement also influences risk. For a midline-extending early T-category tumor, and a clinically N0 ipsilateral neck, the predicted risk is {{< var risk.early_ext_iN0_cN0_II >}}%. This increases to {{< var risk.early_ext_iII+III_cN0_II >}}% when LNLs II and III are involved.  

In summary, midline extension is the primary risk factor for contralateral LNL II involvement, but advanced T-category and extensive ipsilateral involvement also contribute. 

## Contralateral LNL III

As shown in the center panel of @fig-model-risks, the risk of occult contralateral LNL III involvement rarely exceeds 5% and depends strongly on upstream LNL II involvement:  

- For an advanced T-category tumor extending over the midline and with extensive clinical involvement in the ipsilateral LNLs II, III, and IV, but a clinically negative contralateral neck, the risk is only {{< var risk.late_ext_iII+III+IV_cN0_III >}}% (green histogram).  
- If contralateral LNL II is clinically involved, the risk for contralateral LNL III rises to {{< var risk.late_ext_iII_cII_III >}}% (blue histogram).  
- When contralateral LNL II involvement is confirmed by FNA, the risk further increases to {{< var risk.late_ext_iII_cII_FNA_III >}}% (red histogram).  

Even for lateralized tumors, FNA-confirmed involvement in contralateral LNL II predicts a {{< var risk.late_noext_iII_cII_FNA_III >}}% risk for LNL III involvement (orange histogram). This highlights the importance of upstream LNL II involvement in determining the risk for LNL III.

## Contralateral LNL IV

For contralateral LNL IV, the predicted risk is below 1% in most scenarios, even for advanced T-category tumors extending over the midline with extensive ipsilateral and contralateral involvement (green histogram).  

- If contralateral LNL III is also clinically involved, the risk increases to {{< var risk.late_ext_iII+III+IV_cII+III_IV >}}% (blue histogram).  
- When contralateral LNL III involvement is confirmed by FNA, the risk rises significantly to {{< var risk.late_ext_iII+III+IV_cII+III_FNA_IV >}}% (orange histogram).  

This higher risk occurs because FNA confirmation eliminates the possibility of false-positive diagnoses for contralateral LNL III, strongly increasing the likelihood of downstream LNL IV involvement.

## Contralateral LNLs I, V, and VII

Contralateral LNLs I, V, and VII show very low predicted risks for occult disease. Even in advanced T-category tumors extending over the midline with extensive ipsi- and contralateral involvement, the risk for LNLs I and VII remains below 3%.  

For contralateral LNL V, the risk is very small unless there is severe contralateral involvement, including LNL IV, confirmed by FNA. In the extreme case of advanced T-category tumors extending over the mid-sagittal plane, with all ipsilateral LNLs, as well as contralateral LNLs II, III, and IV involved, the risk increases only to {{< var risk.late_ext_iI+II+III+IV+V+VII_cII+III+IV_FNA_V >}}%.

# Discussion {#sec-discussion}

## Summary

This work introduces a formalism to model ipsi- and contralateral lymph node involvement in oropharyngeal SCC patients. Building on a previously developed ipsilateral model [@ludwig_hidden_2021;@ludwig_modelling_2023], we extend it to the contralateral side while preserving the original model's structure. Our extension is both intuitive and interpretable, with parameters learned from a dataset of {{< var data.num_patients >}} patients across four institutions. The model uses clinically diagnosed LNL involvement, tumor T-category, and lateralization to provide personalized risk predictions for occult disease in any LNL of interest.

The model is highly interpretable, with each parameter having a clear, intuitive explanation. Despite its relatively few parameters, it adequately describes the observed data on ipsilateral and contralateral nodal involvement. To our knowledge, this represents the most comprehensive model of lymphatic tumor progression in oropharyngeal SCC, surpassing prior efforts, which were conceptually different, limited in scope, or not trained on real patient data [@benson_markov_2006;@jung_development_2017]. The underlying dataset and code are publicly available, supporting reproducibility and further development.

## Implications for Contralateral Elective Nodal Treatment

Based on a 5% acceptable risk threshold for occult disease in a LNL, the model recommends the following approach for contralateral elective irradiation, assuming the respective LNL is clinically healthy:

- **Lateralized tumors with no contralateral clinical involvement:**  
  Unilateral radiotherapy is sufficient, regardless of T-category or ipsilateral involvement.  
- **Tumors extending over the midline with no contralateral clinical involvement:**  
  Elective irradiation is limited to LNL II.  
- **Contralateral LNL III:**  
  Irradiate LNL III if LNL II is involved, regardless of tumor lateralization, T-category, or ipsilateral involvement. If contralateral LNL II is clinically negative, LNL III is not irradiated unless contralateral LNL IV is involved.  
- **Contralateral LNL IV:**  
  Irradiate only when LNL III involvement is confirmed.  
- **Contralateral LNL V:**  
  Elective irradiation is not recommended in almost all patients. Only in extreme cases, such as advanced T-category tumors with midline extension and confirmed contralateral involvement in LNLs II to IV, irradiation of LNL V may be considered.  
- **Contralateral LNLs I and VII:**  
  Elective irradiation is not recommended unless these levels are clinically involved.

These recommendations align with the model results discussed in @sec-results-risk, providing a basis for refined treatment guidelines. However, they should be interpreted in light of the limitations discussed in @sec-discussion-limitations below.

The model's predictions are already guiding a clinical trial on volume de-escalation at the University Hospital Zurich [@universityofzurich_personalized_2024]. 


## Limitations and Future Work {#sec-discussion-limitations}

### T-Category Dependence

The model uses a single parameter, $p_\text{adv.}$, to account for differences in lymph node involvement patterns between early and advanced T-category tumors. Advanced T-category tumors are modeled as evolving over more time steps, while the probability of spread per time step, governed by the $b_v$ parameters, remains constant. While this approach captures the overall differences between early and advanced T-category tumors well (as seen in @fig-model-prevalences-overall), it is not perfect:

- The observed differences between early and advanced T-category tumors may sometimes be greater or smaller than the model's predictions. This has, for example, been previously noted for ipsilateral LNL I involvement [@ludwig_modelling_2023].
- In the context of the bilateral model, the prevalence of midline extension is overestimated for early T-category tumors and underestimated for advanced T-category tumors, as discussed in @sec-prevalence-midext. 


### Sensitivity and Specificity

As noted in @sec-training-data and further detailed in @sec-consensus, we assumed that the consensus across diagnostic modalities reflects the true state $X_v$ of lymph node involvement. While this assumption is reasonable for pathologically confirmed diagnoses, it is an approximation for clinically diagnosed involvement, which cannot detect occult disease by definition.

The model could, in principle, distinguish between pathologically confirmed and clinically diagnosed involvement by incorporating different sensitivity and specificity values for each diagnostic modality during training. However, this was not implemented in the current work for two reasons:

1. **Simplified evaluation:** This approximation allowed for direct comparison between the observed prevalence of involvement and the model's predictions, enabling an evaluation of the model's ability to describe the data with few interpretable parameters.
2. **Inconsistent literature values:** Reported sensitivity and specificity values for diagnostic modalities show inconsistencies with some observed data, complicating their integration into the model.  

Future efforts could focus on developing methods to rigorously differentiate between pathologically confirmed and clinically diagnosed involvement.

### Tumor Subsites

Oropharyngeal tumors occur in distinct subsites such as the base of the tongue or the tonsils, which may exhibit slightly different lymphatic metastatic spread. Incorporating subsite-specific information into the model could enhance its predictive accuracy. Preliminary investigations suggest that a mixture model may effectively capture these subsite-specific spread patterns [@sarrut_proceedings_2024;@ludwig_modelling_2023a]. 

This approach could also facilitate extending the model to other tumor locations, such as the oral cavity, hypopharynx, and larynx. Including these additional tumor sites would broaden the model's applicability to all HNSCC patients.

# Acknowledgement {.appendix}

We want to thank the research group around Vincent Grégoire from the Centre Léon Bérard in Lyon (France) and Roland Giger's group from the Inselspital in Bern (Switzerland) for collecting, curating, and publishing their data on lymphatic progression patterns. Without their pathologically assessed data on nodal involvement, we could not have trained and validated our model with a large and representative cohort.

Further, this work was supported by:

- the Clinical Research Priority Program "Artificial Intelligence in Oncological Imaging" of the University of Zurich
- the Swiss Cancer Research Foundation under grant number KFS 5645-08-2022

# Consensus on most likely involvement {#sec-consensus .appendix}

The consensus on the most likely involvement state of a LNL was formed as follows: Suppose the involvement status $X_v$ of LNL $v$ was assessed using different diagnostic modalities $\mathcal{O} = \{ \text{MRI}, \text{CT}, \text{pathology}, \ldots \}$, each characterized by their own pair of sensitivity and specificity values $s_N^{\mathcal{o}}$ and $s_P^{\mathcal{o}}$, with $\mathcal{o} \in \mathcal{O}$. These values are tabulated in @tbl-spec-sens. Then we have $|\mathcal{O}|$ observations $z_v^{\mathcal{o}} \in \left[ 0, 1 \right]$, where 0 stands for "healthy" and 1 for "involved". We can then compute the most likely true involvement $X_v$ using the likelihood function

$$
\begin{aligned}
\ell \left( X_v \mid \{ z_v^{\mathcal{o}} \}_{\mathcal{o} \in \mathcal{O}} \right) = \prod_{\mathcal{o} \in \mathcal{O}}
\left( 1 - X_v \right) \cdot &\left[ z_v^{\mathcal{o}} \cdot \left( 1 - s_P^{\mathcal{o}} \right) + \left( 1 - z_v^{\mathcal{o}} \right) \cdot s_P^{\mathcal{o}} \right] \\
+ X_v \cdot &\left[ z_v^{\mathcal{o}} \cdot s_N^{\mathcal{o}} + \left( 1 - z_v^{\mathcal{o}} \right) \cdot (1 - s_N^{\mathcal{o}}) \right]
\end{aligned}
$$

We now assume the true state $X_v$ to take on the value 1 if $\ell \left( X_v = 1 \mid \ldots \right) > \ell \left( X_v = 0 \mid \ldots \right)$ and 0 otherwise. For example, if we have $z_\text{II}^\text{CT} = 0$ and $z_\text{II}^\text{MRI} = 1$ we would compute the following likelihoods:

$$
\begin{aligned}
\ell \left( X_\text{II} = 1 \mid z_\text{II}^\text{CT} = 0, z_\text{II}^\text{MRI} = 1 \right) &= \left( 1 - s_N^\text{CT} \right) \cdot s_N^\text{MRI} = 15.39\% \\
\ell \left( X_\text{II} = 0 \mid z_\text{II}^\text{CT} = 0, z_\text{II}^\text{MRI} = 1 \right) &= s_P^\text{CT} \cdot \left(1 - s_N^\text{MRI}\right) = 14.44\%
\end{aligned}
$$

In this example, we would thus assume the true state to be involved ($X_\text{II} = 1$).

This method of computing a consensus also ensures that the pathology reports always override any conflicting clinical diagnosis, due to pathology's high sensitivity and specificity.

| Modality  | Specificity | Sensitivity |
| :-------- | ----------: | ----------: |
| CT        |         76% |         81% |
| PET       |         86% |         79% |
| MRI       |         63% |         81% |
| FNA       |         98% |         80% |
| pathology |        100% |        100% |

: Specificity and sensitivity values from the literature [@debondt_detection_2007;@kyzas_18ffluorodeoxyglucose_2008]. {#tbl-spec-sens}



# Contralateral Prevalence of Involvement {.appendix}

```{python}
#| echo: false
#| label: tbl-data-strat
#| tbl-cap: Contralateral involvement depending on whether the primary tumor extends over the mid-sagittal plane, the T-category, and how many ipsilateral LNLs were involved.
lnl_cols = get_lnl_cols("contra", lnls=["I", "II", "III", "IV"])
num_ipsi_inv = raw[get_lnl_cols("ipsi")].sum(axis="columns")

contra_inv = raw[lnl_cols].copy()
contra_inv.columns = contra_inv.columns.droplevel([0,1])
contra_inv["t_stage"] = raw[COL.t_stage].apply(lambda x: "early" if x <= 2 else "advanced")
contra_inv["ipsi"] = num_ipsi_inv.map(lambda x: str(x) if x <= 1 else "≥ 2")
contra_inv["midext"] = raw[COL.midext]

grouped = contra_inv.groupby(by=["t_stage", "ipsi", "midext"], dropna=False)
num_involved = grouped.sum()
total = grouped.count()
percent_involved = num_involved / total
idx = total.index.rename(["T-cat.", "ipsi", "Mid. ext."])

total.index = idx
num_involved.index = idx
percent_involved.index = idx

involved = num_involved.join(
  100 * percent_involved,
  rsuffix=" (%)",
)
involved = involved.reindex(sorted(involved.columns), axis="columns")
involved.columns = pd.MultiIndex.from_product(
  [["I", "II", "III", "IV"],
   ["n", "%"]],
  names=["LNL", ""]
)
involved["total", "n"] = total["I"]
(
  involved
  .reset_index()
  .sort_values(
    by=["T-cat.", "ipsi", "Mid. ext."],
    ascending=[False, True, True],
  )
  .style
  .format(precision=2)
  .apply(right_align, subset=involved.columns[0:])
  .apply(highlight_t_stage, subset=["T-cat."])
  .apply(highlight_bool, subset=["Mid. ext."])
  .apply(highlight, subset=["ipsi"], mapping={
    "0": COLORS["green"],
    "1": COLORS["orange"],
    "≥ 2": COLORS["red"],
  })
  .hide()
)
```


# Prevalence of Midline Extension {#sec-prevalence-midext .appendix}

![Comparing the predicted (histograms) and observed (lines depicting beta posteriors) prevalence of midline extension for early (blue) and late (orange) T-category. While the prevalence is predicted correctly when marginalizing over T-category, the model cannot capture the degree of separation observed in the data. Since the tumor's midline extension is virtually always part of the diagnosis and hence *given* when predicting a patient's risk, we do not consider this discrepancy a major issue.](figures/fig_model_prevalences_midext.svg){#fig-model-prevalences-midext width="50%"}

In @fig-model-prevalences-midext, we plot the prevalence of midline extension in the data versus our model's prediction. The model adjusts the parameter $p_\epsilon$ to correctly predict the overall proportion of patients with midline extension, but it cannot match the large spread between early and advanced T-category seen in the data. To achieve that, the model would need to increase $p_\text{adv.}$ and decrease $p_\epsilon$. But since the parameter $p_\text{adv.}$ also determines the differences in LNL involvement between early and advanced T-category, the model does not have that freedom.

However, we do not consider this discrepancy a major limitation of the model: In a practical application of the model we are not interested in the probability of midline extension, as it is always possible to assess it with high certainty for the patient at hand. That is also the reason why we initially modelled the midline extension *not* as a random variable, but as a global risk factor that would have been turned on or off from the onset of a patient's disease evolution. This, however, lead to overly high risks for contralateral involvement in advanced T-category patients with midline extension, because then the model assumes an increased spread to the contralateral side from the onset of the disease. Which is probably not true in a majority of those cases. Thus, treating it as a random variable that only becomes true during a patient's disease evolution resulted in a better description of the data.

Formally, the wrong prediction of midline extension prevalence makes little difference, since it is always given: Instead of $P\left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \mid \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \right)$, we typically compute $P\left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \mid \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c}, \epsilon \right)$, which does not suffer from the wrong probability of midline extension, as the distribution over hidden states is renormalized:

$$
P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c} \mid \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c}, \epsilon \right) = \frac{P \left( \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c} \mid \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \right) P \left( \mathbf{X}^\text{i}, \mathbf{X}^\text{c}, \epsilon \right)}{P \left( \mathbf{Z}^\text{i}, \mathbf{Z}^\text{c}, \epsilon \right)}
$$

Note that a distribution over $\epsilon$ appears both in the enumerator and the denominator, which largely cancel each other, leaving only the midline extension's effect on the distribution over hidden states in the prediction.
