Probabilistic Model of Bilateral Lymphatic Spread in Head and Neck Cancer

In [2]:

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)

In [3]:

# 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) (Grégoire et al. 2003, 2014, 2018; Eisbruch et al. 2002; Biau et al. 2019; Chao et al. 2002; Vorwerk and Hess 2011; Ferlito, Silver, and Rinaldo 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) (Grégoire et al. 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 et al. 2019).

To better quantify individualized risk of occult disease, we previously developed an intuitive probabilistic hidden Markov model (HMM) (Ludwig et al. 2021, 2023), originally based on a conceptually similar a Bayesian network model (Pouymayou et al. 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. Section 2 presents a multi-centric dataset on lymph node involvement in 833 OPSCC patients, identifying key risk factors for contralateral lymph node involvement and outlining requirements for a bilateral model extension (section 2.4).

2. Section 4 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 section 5.

3. Section 6 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 section 8.

Data on Lymphatic Progression Patterns

To develop models for lymphatic tumor progression for all relevant LNLs, including contralateral regions, we compiled a detailed dataset of 833 patients with newly diagnosed oropharyngeal squamous cell carcinomas (OPSCC) (Ludwig et al. 2022; Ludwig, Schubert, Barbatei, Bauwens, Werlen, et al. 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 table 1.

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.

In [4]:

In [4]:

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

Table 1: 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.

Institution	Total	Age (median)	Neck Dissection	N0	Early T-Cat.	Mid. Ext.
Centre Léon Bérard	325	60	100%	19%	69%	18%
Inselspital Bern	74	61	100%	18%	66%	14%
University Hospital Zurich	287	66	26%	18%	52%	31%
Vall d'Hebron Barcelona Hospital	147	58	5%	21%	34%	34%

Consensus on Involvement Status

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 (De Bondt et al. 2007; Kyzas et al. 2008).

The consensus process is detailed in section 10. 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, where it is also available for interactive exploration, from GitHub, or from Zenodo. Additional details about the datasets and data format are provided in designated Data-in-Brief publications (Ludwig et al. 2022; Ludwig, Schubert, Barbatei, Bauwens, Werlen, et al. 2024). These publications do not include the latest data from HVH, which will be described in a future publication.

Patterns of Contralateral Involvement

The datasets enable analysis of correlations between contralateral LNL involvement and key risk factors. In figure 1, 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.

Figure 1: 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.

Midline Extension

The bottom left panel of figure 1 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.

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.

Correlation of Risk Factors

Midline extension, T-category, and ipsilateral involvement are interrelated risk factors for contralateral metastasis. For instance, 45.6% of advanced T-category tumors exhibit midline extension compared to 6.1% 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 figure 1 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 1.2% (1 of 86 patients) show involvement in contralateral level II. This proportion increases to 8.8% (24 of 272) if ipsilateral level II is involved, to 15.7% (14 of 89) if ipsilateral levels II and III are involved, and further to 22.2% (12 of 54) for advanced T-category tumors with ipsilateral levels II and III involved.

Requirements for a Bilateral Model

Based on the observations in section 2.3 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 et al. 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

This paper builds on the previously developed unilateral model for ipsilateral lymph node involvement presented in (Ludwig, Schubert, Barbatei, Bauwens, Hoffmann, et al. 2024). Below, we briefly recap the unilateral model to introduce the notation required for extending the framework to a bilateral model, described in section 4. For further details on the ipsilateral model, we refer to earlier publications (Ludwig et al. 2021; Ludwig, Schubert, Barbatei, Bauwens, Hoffmann, et al. 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\} \tag{1}\]

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) \tag{2}\] 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) \tag{3}\] 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 figure 2 which was obtained by maximizing the model evidence as described in (Ludwig, Schubert, Barbatei, Bauwens, Hoffmann, et al. 2024).

Figure 2: 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, …).

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} \tag{4}\]

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) \tag{5}\]

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 equation 5 depends on the spread parameters shown in figure 2 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)} \tag{6}\]

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 equation 5. 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 \tag{7}\]

The methodology we use to infer the model’s parameters is detailed in section 5.2.

Extension to a Bilateral Model

A straightforward approach to modeling contralateral lymphatic spread would be to use two independent unilateral models, as described in section 3, 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 section 2.3, particularly the observed increase in contralateral involvement with greater severity of ipsilateral spread.

Thus, we extend the formalism in section 3 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)} \tag{8}\]

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} \tag{9}\]

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 equation 9 along with equation 4, 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} \tag{10}\]

Parameter Symmetries

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 figure 2.
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} \tag{11}\] for all \(v \leq V\) and \(r \in \operatorname{pa}(v)\) being all nodes that spread to \(v\).
3. 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}} \tag{12}\] 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)} \tag{13}\]

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

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 equation 4, 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 equation 12. 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 equation 4:

\[ \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 equation 10, 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} \tag{14}\]

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)} \tag{15}\]

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}} \tag{16}\]

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 equation 3. 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} \tag{17}\]

Computational Methods

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

Training Data

We trained the model using the dataset of 833 patients described in section 2. 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

We used the Python package emcee (Foreman-Mackey et al. 2013) for parameter inference, implementing efficient MCMC sampling with parallel affine-invariant samplers. The sampling algorithms employed differential evolution moves (Ter Braak and Vrugt 2008; Nelson, Ford, and Payne 2013), with the likelihood implemented by our lymph-model Python package.

We initialized 12 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 5.0e-2.
2. The autocorrelation estimate dropped below \(n\) / 50, 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 10 additional samples, spaced 10 steps apart.

We verified sampling convergence in figure 3 by examining the MCMC chain’s autocorrelation time and walker acceptance fractions.

Figure 3: 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.

Comparing the Observed and Predicted Prevalence of Involvement Patterns

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 figure 7 (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=\) 379 such patients were observed, with \(k=\) 27 exhibiting contralateral LNL II involvement – an observed prevalence of \(q=\) 7.1%. 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 figure 7, we consider the case of early T-category tumors extending over the midline. The dataset contains only 29 such patients, out of which 6 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 equation 17, 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

In table 2, 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 et al. 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\) (33.9%) 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}\) (14.2%) compared to the value of \(b_3^\text{i}\) (5.5%) shows that the ipsilateral LNL III is rarely involved without the upstream involvement of LNL II.

In [5]:

In [5]:

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

Table 2: 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\).

Parameter	Mean	Std. Dev.
Mid. ext. probability	8.16%	± 0.48%
ipsi: T ➜ I	2.80%	± 0.26%
ipsi: T ➜ II	34.89%	± 1.40%
ipsi: T ➜ III	5.45%	± 0.66%
ipsi: T ➜ IV	0.94%	± 0.18%
ipsi: T ➜ V	1.83%	± 0.22%
ipsi: T ➜ VII	2.32%	± 0.26%
contra: T ➜ I	0.29%	± 0.09%
contra: T ➜ II	2.46%	± 0.29%
contra: T ➜ III	0.14%	± 0.07%
contra: T ➜ IV	0.19%	± 0.08%
contra: T ➜ V	0.05%	± 0.04%
contra: T ➜ VII	0.50%	± 0.17%
Mixing ⍺	33.87%	± 4.32%
I ➜ II	62.50%	± 16.69%
II ➜ III	14.23%	± 1.64%
III ➜ IV	15.86%	± 1.93%
IV ➜ V	14.58%	± 3.76%
late T-cat. binom. prob.	44.98%	± 1.99%

Illustration of the model

In this subsection, we illustrate key aspects of the mathematical framework introduced earlier. The top panel of figure 4 shows the prior distribution over diagnosis times, \(P(t)\). Based on the parameterization, early T-category tumors are on average diagnosed after 3 time steps, while advanced T-category tumors are diagnosed later, averaging 4.5 time steps. This is due to the learned value of \(p_\text{adv.}\), which is 45% for advanced T-category tumors. The tumor’s average probability per time step of growing over the midline, \(p_\epsilon\), was found to be 8.2%. 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 figure 4). 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.

Figure 4: 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.

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 figure 5, which represents the calculation defined in equation 14. 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 table 2.

Figure 5: Visual representation of equation 14. 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 figure 6.

The left matrix in figure 5 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 figure 6.

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.

Figure 6: 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 equation 14, visualized in figure 5.

Prevalence predictions for contralateral involvement

The bilateral model was designed to meet the requirements outlined in section 2.4. 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 figure 7, we compare the prevalence of contralateral involvement for LNLs II, III, and IV.

Figure 7: 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.

Figure 7 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 7.1% for early T-category lateralized tumors to 39.2% for advanced T-category tumors that cross the midline. Similarly, contralateral LNL III involvement rises from around 1.5% for early T-category lateralized tumors to nearly 14.2% for advanced T-category midline-extending tumors.

Influence of Upstream Involvement on Contralateral Metastasis

Figure 8: 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).

Figure 8 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 164 patients with an advanced tumor crossing the midline and involvement in the upstream LNL II, 19 had involvement in the contralateral LNL III. In contrast, out of 379 patients with a clearly lateralized early T-category tumor and no upstream involvement, only 1 showed contralateral LNL III involvement. This is well captured by our model.

Correlation between Ipsi- and Contralateral Involvement

Figure 9: 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).

In figure 9, 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

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% (De Bondt et al. 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.

Figure 10: 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%).

Contralateral LNL II

The predicted risk of occult disease in contralateral LNL II is shown in the left panel of figure 10. 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 1.6% 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 7.6% (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 9.1% risk of contralateral LNL II disease (red histogram).
  
• For the same scenario but an advanced T-category tumor, the risk rises to 11.3% (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 7.6%. This increases to 9.1% 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 figure 10, 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 2.1% (green histogram).
  
• If contralateral LNL II is clinically involved, the risk for contralateral LNL III rises to 4.8% (blue histogram).
  
• When contralateral LNL II involvement is confirmed by FNA, the risk further increases to 6.9% (red histogram).

Even for lateralized tumors, FNA-confirmed involvement in contralateral LNL II predicts a 5.6% 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 2.8% (blue histogram).
  
• When contralateral LNL III involvement is confirmed by FNA, the risk rises significantly to 5.7% (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 5.5%.

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 et al. 2021, 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 833 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, Whipple, and Kalet 2006; Jung et al. 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 section 7, providing a basis for refined treatment guidelines. However, they should be interpreted in light of the limitations discussed in section 8.3 below.

The model’s predictions are already guiding a clinical trial on volume de-escalation at the University Hospital Zurich (University of Zurich 2024).

Limitations and Future Work

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 figure 7), 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 et al. 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 section 12.

Sensitivity and Specificity

As noted in section 5.1 and further detailed in section 10, 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 and Rit 2024; Ludwig 2023).

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

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

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 table 3. 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.

Table 3: Specificity and sensitivity values from the literature (De Bondt et al. 2007; Kyzas et al. 2008).

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

Contralateral Prevalence of Involvement

In [6]:

In [6]:

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

Table 4: Contralateral involvement depending on whether the primary tumor extends over the mid-sagittal plane, the T-category, and how many ipsilateral LNLs were involved.

T-cat.	ipsi	Mid. ext.	I		II		III		IV		total
			n	%	n	%	n	%	n	%	n
early	0	False	0	0.00	1	1.16	0	0.00	0	0.00	86
early	0	True	0	0.00	1	10.00	1	10.00	0	0.00	10
early	0	nan	0	0.00	0	0.00	0	0.00	0	0.00	12
early	1	False	1	0.53	11	5.82	2	1.06	1	0.53	189
early	1	True	1	11.11	2	22.22	0	0.00	0	0.00	9
early	1	nan	0	0.00	3	9.68	0	0.00	1	3.23	31
early	≥ 2	False	1	0.96	15	14.42	3	2.88	3	2.88	104
early	≥ 2	True	0	0.00	3	30.00	4	40.00	1	10.00	10
early	≥ 2	nan	0	0.00	3	13.64	0	0.00	0	0.00	22
advanced	0	False	0	0.00	2	5.88	0	0.00	0	0.00	34
advanced	0	True	0	0.00	3	12.50	0	0.00	0	0.00	24
advanced	0	nan	0	0.00	0	0.00	0	0.00	0	0.00	3
advanced	1	False	0	0.00	3	4.55	0	0.00	0	0.00	66
advanced	1	True	1	1.64	18	29.51	5	8.20	1	1.64	61
advanced	1	nan	0	0.00	0	0.00	0	0.00	0	0.00	9
advanced	≥ 2	False	4	5.26	16	21.05	7	9.21	4	5.26	76
advanced	≥ 2	True	3	3.80	46	58.23	18	22.78	8	10.13	79
advanced	≥ 2	nan	0	0.00	0	0.00	0	0.00	0	0.00	8

Prevalence of Midline Extension

Figure 11: 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.

In figure 11, 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.

Benson, Noah, Mark Whipple, and Ira J Kalet. 2006. “A Markov Model Approach to Predicting Regional Tumor Spread in the Lymphatic System of the Head and Neck.” AMIA ... Annual Symposium Proceedings. AMIA Symposium, 31–35.

Biau, Julian, Michel Lapeyre, Idriss Troussier, Wilfried Budach, Jordi Giralt, Cai Grau, Joanna Kazmierska, et al. 2019. “Selection of Lymph Node Target Volumes for Definitive Head and Neck Radiation Therapy: A 2019 Update.” Radiotherapy and Oncology 134 (May): 1–9. https://doi.org/10.1016/j.radonc.2019.01.018.

Chao, K. S.Clifford, Franz J Wippold, Gokhan Ozyigit, Binh N Tran, and James F Dempsey. 2002. “Determination and Delineation of Nodal Target Volumes for Head-and-Neck Cancer Based on Patterns of Failure in Patients Receiving Definitive and Postoperative IMRT.” International Journal of Radiation Oncology*Biology*Physics 53 (5): 1174–84. https://doi.org/10.1016/S0360-3016(02)02881-X.

De Bondt, R. B. J., P. J. Nelemans, P. A. M. Hofman, J. W. Casselman, B. Kremer, J. M. A. Van Engelshoven, and R. G. H. Beets-Tan. 2007. “Detection of Lymph Node Metastases in Head and Neck Cancer: A Meta-Analysis Comparing US, USgFNAC, CT and MR Imaging.” European Journal of Radiology 64 (2): 266–72. https://doi.org/10.1016/j.ejrad.2007.02.037.

Eisbruch, Avraham, Robert L. Foote, Brian O’Sullivan, Jonathan J. Beitler, and Bhadrasain Vikram. 2002. “Intensity-Modulated Radiation Therapy for Head and Neck Cancer: Emphasis on the Selection and Delineation of the Targets.” Seminars in Radiation Oncology 12 (3): 238–49. https://doi.org/10.1053/srao.2002.32435.

Ferlito, Alfio, Carl E. Silver, and Alessandra Rinaldo. 2009. “Elective Management of the Neck in Oral Cavity Squamous Carcinoma: Current Concepts Supported by Prospective Studies.” British Journal of Oral and Maxillofacial Surgery 47 (1): 5–9. https://doi.org/10.1016/j.bjoms.2008.06.001.

Foreman-Mackey, Daniel, David W. Hogg, Dustin Lang, and Jonathan Goodman. 2013. “Emcee: The MCMC Hammer.” Publications of the Astronomical Society of the Pacific 125 (925): 306–12. https://doi.org/10.1086/670067.

Grégoire, Vincent, Kian Ang, Wilfried Budach, Cai Grau, Marc Hamoir, Johannes A. Langendijk, Anne Lee, et al. 2014. “Delineation of the Neck Node Levels for Head and Neck Tumors: A 2013 Update. DAHANCA, EORTC, HKNPCSG, NCIC CTG, NCRI, RTOG, TROG Consensus Guidelines.” Radiotherapy and Oncology 110 (1): 172–81. https://doi.org/10.1016/j.radonc.2013.10.010.

Grégoire, Vincent, Mererid Evans, Quynh-Thu Le, Jean Bourhis, Volker Budach, Amy Chen, Abraham Eisbruch, et al. 2018. “Delineation of the Primary Tumour Clinical Target Volumes (CTV-P) in Laryngeal, Hypopharyngeal, Oropharyngeal and Oral Cavity Squamous Cell Carcinoma: AIRO, CACA, DAHANCA, EORTC, GEORCC, GORTEC, HKNPCSG, HNCIG, IAG-KHT, LPRHHT, NCIC CTG, NCRI, NRG Oncology, PHNS, SBRT, SOMERA, SRO, SSHNO, TROG Consensus Guidelines.” Radiotherapy and Oncology 126 (1): 3–24. https://doi.org/10.1016/j.radonc.2017.10.016.

Grégoire, Vincent, Peter Levendag, Kian K. Ang, Jacques Bernier, Marijel Braaksma, Volker Budach, Cliff Chao, et al. 2003. “CT-based Delineation of Lymph Node Levels and Related CTVs in the Node-Negative Neck: DAHANCA, EORTC, GORTEC, NCIC,RTOG Consensus Guidelines.” Radiotherapy and Oncology 69 (3): 227–36. https://doi.org/10.1016/j.radonc.2003.09.011.

Jung, Hyunggu, Anthony Law, Eli Grunblatt, Lucy L. Wang, Aaron Kusano, Jose L. V. Mejino, and Mark E. Whipple. 2017. “Development of a Novel Markov Chain Model for the Prediction of Head and Neck Squamous Cell Carcinoma Dissemination.” AMIA Annual Symposium Proceedings 2016 (February): 1832–39.

Kyzas, Panayiotis A., Evangelos Evangelou, Despina Denaxa-Kyza, and John P. A. Ioannidis. 2008. “18F-fluorodeoxyglucose Positron Emission Tomography to Evaluate Cervical Node Metastases in Patients with Head and Neck Squamous Cell Carcinoma: A Meta-Analysis.” JNCI: Journal of the National Cancer Institute 100 (10): 712–20. https://doi.org/10.1093/jnci/djn125.

Ludwig, Roman. 2023. “Modelling Lymphatic Metastatic Progression in Head and Neck Cancer.” Ludwig, Roman. Modelling Lymphatic Metastatic Progression in Head and Neck Cancer. 2023, University of Zurich, Faculty of Science. PhD thesis, Zurich: University of Zurich. https://doi.org/10.5167/uzh-231470.

Ludwig, Roman, Jean-Marc Hoffmann, Bertrand Pouymayou, Grégoire Morand, Martina Broglie Däppen, Matthias Guckenberger, Vincent Grégoire, Panagiotis Balermpas, and Jan Unkelbach. 2022. “A Dataset on Patient-Individual Lymph Node Involvement in Oropharyngeal Squamous Cell Carcinoma.” Data in Brief 43 (August): 108345. https://doi.org/10.1016/j.dib.2022.108345.

Ludwig, Roman, Bertrand Pouymayou, Panagiotis Balermpas, and Jan Unkelbach. 2021. “A Hidden Markov Model for Lymphatic Tumor Progression in the Head and Neck.” Scientific Reports 11 (1): 12261. https://doi.org/10.1038/s41598-021-91544-1.

Ludwig, Roman, Adrian Daniel Schubert, Dorothea Barbatei, Lauence Bauwens, Jean-Marc Hoffmann, Sandrine Werlen, Olgun Elicin, et al. 2024. “Modelling the Lymphatic Metastatic Progression Pathways of OPSCC from Multi-Institutional Datasets.” Scientific Reports 14 (1): 15750. https://doi.org/10.1038/s41598-024-66012-1.

Ludwig, Roman, Adrian Schubert, Dorothea Barbatei, Lauence Bauwens, Jean-Marc Hoffmann, Sandrine Werlen, Olgun Elicin, et al. 2023. “Modelling the Lymphatic Metastatic Progression Pathways of OPSCC from Multi-Institutional Datasets.” arXiv. https://doi.org/10.48550/arXiv.2312.11270.

Ludwig, Roman, Adrian Schubert, Dorothea Barbatei, Laurence Bauwens, Sandrine Werlen, Olgun Elicin, Matthias Dettmer, et al. 2024. “A Multi-Centric Dataset on Patient-Individual Pathological Lymph Node Involvement in Head and Neck Squamous Cell Carcinoma.” Data in Brief 52 (February): 110020. https://doi.org/10.1016/j.dib.2023.110020.

Nelson, Benjamin, Eric B. Ford, and Matthew J. Payne. 2013. “RUN DMC: AN EFFICIENT, PARALLEL CODE FOR ANALYZING RADIAL VELOCITY OBSERVATIONS USING N -BODY INTEGRATIONS AND DIFFERENTIAL EVOLUTION MARKOV CHAIN MONTE CARLO.” The Astrophysical Journal Supplement Series 210 (1): 11. https://doi.org/10.1088/0067-0049/210/1/11.

Pouymayou, Bertrand, Panagiotis Balermpas, Oliver Riesterer, Matthias Guckenberger, and Jan Unkelbach. 2019. “A Bayesian Network Model of Lymphatic Tumor Progression for Personalized Elective CTV Definition in Head and Neck Cancers.” Physics in Medicine & Biology 64 (16): 165003. https://doi.org/10.1088/1361-6560/ab2a18.

Sarrut, David, and Simon Rit. 2024. Proceedings of the XX-th International Conference on the Use of Computers in Radiation Therapy (ICCR). XX-Th International Conference on the Use of Computers in Radiation Therapy (ICCR). Lyon, France.

Ter Braak, Cajo J. F., and Jasper A. Vrugt. 2008. “Differential Evolution Markov Chain with Snooker Updater and Fewer Chains.” Statistics and Computing 18 (4): 435–46. https://doi.org/10.1007/s11222-008-9104-9.

University of Zurich. 2024. “Personalized Volume-deescalated Elective Nodal Irradiation in Oropharyngeal Head and Neck.” Clinical Trial Registration NCT06563362. clinicaltrials.gov.

Vorwerk, Hilke, and Clemens F Hess. 2011. “Guidelines for Delineation of Lymphatic Clinical Target Volumes for High Conformal Radiotherapy: Head and Neck Region.” Radiation Oncology 6 (1): 97. https://doi.org/10.1186/1748-717X-6-97.