# Tumor kinetics modeling outperforms best overall response for prediction of survival in head and neck squamous cell carcinoma

Kevin Atsou, Anne Auperin, Joël Guigay, Sébastien Salas, Sébastien Benzekry

► To cite this version:

Kevin Atsou, Anne Auperin, Joël Guigay, Sébastien Salas, Sébastien Benzekry. Tumor kinetics modeling outperforms best overall response for prediction of survival in head and neck squamous cell carcinoma. 2024. hal-04558029

HAL Id: hal-04558029

https://inria.hal.science/hal-04558029

Preprint submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mechanistic Learning for Predicting Survival Outcomes in Head and Neck Squamous Cell Carcinoma

Kevin Atsou, Anne Auperin, Jôel Guigay, Sébastien Salas, Sébastien Benzekry

► To cite this version:

Kevin Atsou, Anne Auperin, Jôel Guigay, Sébastien Salas, Sébastien Benzekry. Mechanistic Learning for Predicting Survival Outcomes in Head and Neck Squamous Cell Carcinoma. CPT: Pharmacometrics and Systems Pharmacology, 2024, 10.1002/psp4.13294. hal-04558029v2

HAL Id: hal-04558029

https://inria.hal.science/hal-04558029v2

Submitted on 29 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tumor kinetics modeling outperforms best overall response for prediction of survival in head and neck squamous cell carcinoma

Kevin Atsou¹, Anne Auperin², Jôel Guigay³, Sébastien Salas¹,⁴, Sebastien Benzekry ¹, *

1. COMPUTational pharmacology and clinical Oncology Department, Inria Sophia Antipolis – Méditerranée, Cancer Research Center of Marseille, Inserm UMR1068, CNRS UMR7258, Aix Marseille University UM105, Marseille, France.

2. Biostatistical and Epidemiological Division, Institut Gustave Roussy, Villejuif, France.

3. Clinical Oncology, Centre Antoine Lacassagne, Nice, France.

4. Clinical Oncology, Hôpital Timone, Aix-Marseille University, Marseille, France.

The authors declare no conflict of interest.

Social media handles:

Authors: @SBenzekry, @KevinAtsou

Institutions : @aphm_ actu, @inria_sophia, @crcm_marseille
ARTICLE OPEN ACCESS

# Mechanistic Learning for Predicting Survival Outcomes in Head and Neck Squamous Cell Carcinoma

Kevin Atsou¹ | Anne Auperin² | Jôel Guigay³ | Sébastien Salas¹,⁴ | Sebastien Benzekry¹

¹COMPutational Pharmacology and Clinical Oncology Department, Inria Sophia Antipolis - Méditerranée, Cancer Research Center of Marseille, Inserm UMR1068, CNRS UMR7258, Aix Marseille University UM105, Marseille, France | ²Biostatistical and Epidemiological Division, Institut Gustave Roussy, Villejuif, France | ³Clinical Oncology, Centre Antoine Lacassagne, Nice, France | ⁴Clinical Oncology, Hôpital Timone, Aix-Marseille University, Marseille, France

**Correspondence:** Kevin Atsou (atsoukevin@gmail.com; kokoukevin.atsou@pfizer.com) | Sebastien Benzekry (sebastien.benzekry@inria.fr)

**Received:** 7 June 2024 | **Revised:** 18 November 2024 | **Accepted:** 5 December 2024

**Funding:** This work is part of the QUANTIC project funded by ITMO Cancer AVIESAN and the French Institute National du Cancer (grant #19CM148-00) and of the DIGPHAT project which was supported by a grant by the French government, managed by the National Research Agency (ANR), under the France 2030 program, reference ANR-22-PESN-0017.

**Keywords:** deep learning | head and neck squamous cell carcinoma | machine learning | survival analysis | tumor kinetics

## ABSTRACT

We employed a mechanistic learning approach, integrating on-treatment tumor kinetics (TK) modeling with various machine learning (ML) models to address the challenge of predicting post-progression survival (PPS)—the duration from the time of documented disease progression to death—and overall survival (OS) in Head and Neck Squamous Cell Carcinoma (HNSCC). We compared the predictive power of model-derived TK parameters versus RECIST and assessed the efficacy of nine TK-OS ML models against conventional survival models. Data from 526 advanced HNSCC patients treated with chemotherapy and cetuximab in the TPExtreme trial were analyzed using a double-exponential model. TK parameters from the first line and maintenance (TKL1) or after four cycles (TK4) were used to predict PPS and post-cycle 4 OS (OS4), combined with 12 baseline parameters. While ML algorithms underperformed compared to the Cox model for PPS, a random survival forest was superior for OS prediction using TK4 and surpassed RECIST-based metrics. This model demonstrated unbiased OS4 prediction, suggesting its potential for improving HNSCC treatment evaluation.

**Trial Registration:** ClinicalTrials.gov identifier: NCT02268695.

## 1 | Introduction

Head and Neck Squamous Cell Carcinomas (HNSCCs) arise from the mucosal epithelium of the oral cavity, pharynx and the larynx, posing a significant clinical challenge due to their aggressive nature and associated morbidity. The standard of care for resectable HNSCCs involves surgery followed by radiotherapy, with or without chemotherapy [1, 2]. For unresectable loco-regional HNSCCs, chemo-radiotherapy is the preferred treatment. In cases of local or metastatic recurrence where patients are not eligible for local treatments

such as surgery or radiotherapy, standard treatment options includes the EXTREME regimen, which combines cisplatin and 5-fluorouracil (5-FU) with cetuximab [3] or a combination of cisplatin/carboplatin and 5-fluorouracil (5-FU) with pembrolizumab depending on the PD-L1 status [4]. These regimens are typically administered over 6 cycles of 21 days each, as a first line therapeutic approach. The GORTEC 2014-01 TPExtreme clinical trial aimed to compare the efficacy and safety of the TPEx regimen—4 cycles of docetaxel in combination with cisplatin and cetuximab followed by 500 mg/m² Cetuximab maintenance every 2 weeks—with the EXTREME

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

© 2024 The Author(s). CPT: Pharmacometrics & Systems Pharmacology published by Wiley Periodicals LLC on behalf of American Society for Clinical Pharmacology and Therapeutics.
# ABSTRACT

## Background

In the realm of Head and Neck Squamous Cell Carcinoma (HNSCC) treatment, accurately predicting post-progression survival (PPS) and overall survival (OS) from response metrics remains a pivotal challenge. The current mainstay for response evaluation is based on the response evaluation criteria in solid tumors (RECIST) that do not account for the entire tumor response kinetics (TK). Model-derived TK parameters have been reported to carry significant OS predictive power but a rigorous comparison of such metrics versus RECIST remains to be established, and TK-OS often suffer from time-dependent covariate bias. Moreover, the comparative efficacy of TK-OS machine learning (ML) versus conventional survival models remains to be assessed.

## Methods

Data on 526 advanced HNSCC patients treated with chemotherapy and cetuximab were collected from the TPExtreme clinical trial. A double-exponential model was selected to describe both first line and maintenance treatment periods. Model-based TK parameters derived from data available after first line and maintenance (TKL1) or after 4 treatment cycles (TK4) were considered to respectively predict post-progression survival (PPS) and post-4 cycles survival (OS4). They were combined with 13 baseline (BSL) clinical parameters and integrated into nine survival ML algorithms that were benchmarked for their predictive performances. A training set and a test set with 70/30 proportions were defined for model calibration and evaluation, respectively.

## Results

Using TKL1 for PPS prediction, the ML algorithms felt short the classical Cox proportional hazards model. For OS4 prediction using TK4, a random survival forest (RSF) emerged as the best model. Besides the performance status that emerged as the most important, the next four top features were TK4 metrics. Furthermore, TK4 metrics outperformed RECIST-based best overall response (BOR) (c-index of 0.55 (0.52 - 0.57) vs 0.58 (0.57 - 0.6) on the training set, *p* < 0.0001 and 0.5 (0.47 - 0.52) vs 0.53 (0.50 - 0.55) on the test set, *p* < 0.0001, BOR + BSL vs TK4 + BSL, respectively).

## Conclusion

The integration of TK4 model-based parameters with a random survival algorithm provides unbiased and superior OS4 predictions than RECIST in HNSCC treated with chemotherapy.
## Summary

* What is the current knowledge on the topic?

* Current surrogate markers for predicting overall survival (OS) in head and neck squamous cell carcinoma (HNSCC) rely on RECIST criteria. However, these do not account for the full dynamics of tumor kinetics (TK), limiting their predictive accuracy.

* What question did this study address?

* This study evaluated the efficacy of model-derived TK parameters versus RECIST-based metrics for predicting post-progression survival (PPS) and overall survival (OS) in HNSCC patients, comparing traditional survival models with advanced machine learning (ML) and deep learning (DL) approaches.

* What does this study add to our knowledge?

* The study demonstrated that TK parameters, especially when combined with random survival forests (RSF), provided superior OS predictions compared to RECIST metrics. While complex ML and DL models did not significantly outperform classical models, RSF was the most effective for OS prediction using TK data.

* How might this change drug discovery, development, and/or therapeutics?

* Integrating TK-based metrics with ML models could improve early predictions of individual patient outcomes, potentially leading to more personalized and effective treatment strategies in HNSCC. This approach could refine the evaluation and monitoring of treatment responses beyond traditional RECIST guidelines.

regimen, potentially offering an alternative for the standard of care in the first-line treatment, particularly in patients who may not be good candidates for up-front pembrolizumab treatment [5]. Despite these therapeutic strides, the grim reality of a median survival oscillating between 12 and 15 months post-treatment [5] underlines a pressing need for more precise prognostic models to guide clinical decision-making and personalize therapeutic strategies.

Accurate prediction of overall survival (OS) from early surrogate markers is crucial for both for clinical care and drug development, especially to forecast the outcomes of phase 3 trials based on early interim data. Current surrogate markers to predict OS are based on the response evaluation criteria in solid tumors (RECIST) [6]. RECIST defines the best overall response (BOR) as the most favorable response recorded from treatment initiation to disease progression, based on the relative change in the sum of the largest diameters (SLD) of target lesions from baseline. Additionally, RECIST allows the definition of progression-free survival (PFS), which is the time from the start of treatment until objective tumor progression or death. However, PFS and BOR have been weakly supported as surrogate markers for OS [7–11].

To enhance prognostic models, on-treatment longitudinal SLD measurements have been employed to model TK using semi-mechanistic “tumor growth inhibition” (TGI) models [12, 13]. These models provide insights into the dynamics of tumor response

to treatment over time and generate TK parameters for prognostic modeling. Several studies combining estimated TK parameters to survival analysis (TK-OS) revealed the importance of TK parameters in predicting OS, for example, in metastatic colorectal cancer or non-small cell lung cancer (NSCLC) [12–14]. However, a direct comparison between the predictive abilities of TK parameters versus RECIST-based BOR has not yet been performed.

Initially introduced to assist drug development (e.g., predicting phase 3 outcomes from phase 2 data) [12], TK parameters can also be used for individual OS prediction. However, in such context, on-treatment longitudinal data cannot be used to predict survival at a pre-treatment time point, as is the case when the outcome is OS [15, 16]. This phenomenon has been referred to as a type of time-dependent covariate bias and used as a criticism of some TK-OS studies [16]. For instance, using the entire duration of on-treatment data to predict OS is problematic because the total on-treatment duration is itself a strong predictor of OS. One possibility to alleviate this bias is to identify TK parameters from SLD measurements within a restricted timeframe to predict subsequent survival. Predicting post-progression survival (PPS), defined as the time from disease progression until death or last follow-up, is also a suitable alternative to OS. Specifically, we focused on two well-defined predictive tasks: (1) predicting post-cycle 4 OS (OS4) from tumor kinetics data up to that point (TK4-OS4) and (2) predicting post-progression survival (PPS) from the first line tumor kinetics only (TK1L-PPS). Note that point (2) is of particular interest to clinical oncologists because the time of failure of a first treatment is pivotal for therapeutic decision. Being able to accurately predict PPS could help to choose whether to give a treatment or not, and what treatment to give as second line in the first case.

To address this, we harnessed the strengths of population TK modeling, utilizing non-linear mixed effects (NLME), coupled with machine learning (ML) survival analysis. By benchmarking the currently available ML survival models, we compared the predictive power of TK parameters against BOR for survival prediction.

## 2 | Methods

### 2.1 | Patients

Data consisted of all patients enrolled in the GORTEC 2014-01 TPExtreme clinical trial [5]. Included patients had at least one measurable lesion according to the Response Evaluation Criteria in Solid Tumors (RECIST) v1.1 [6], ECOG performance status of 1 or 0, and were aged 18–70 years with histologically confirmed metastatic or recurrent squamous-cell carcinoma of the oral cavity, hypopharynx, oropharynx, or larynx, unsuitable for locoregional curative treatments. A total of 526 participants were randomly assigned (1:1) to each treatment group. On one hand, the TPEx regimen was administrated and consisted of 4 cycles of a combination of docetaxel, cisplatin, and cetuximab. On the other hand, the EXTREME regimen consisted of six cycles of the standard of care (fluorouracil, cisplatin, and cetuximab) (Figure 1A).

All enrolled patients gave written informed consent before any study procedure. The trial was carried out in respect of good clinical practice guidelines and the declaration of Helsinki. It
**Keywords:** Survival analysis, Machine learning, Tumor kinetics, Head and neck squamous cell carcinoma
A

<figure><img src="image_2.png" /></figure>

B

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Treatment history and kinetics-ML prediction of survival. (A) Treatment history of a sample of patients. Cetuximab was the treatment given during maintenance. (B) Steps describing mechanistic learning. It combines modeling tumor kinetics and ML for survival prediction.</figcaption></figure>

was approved by competent authorities and ethics committees
in France, Spain, and Germany and registered with [ClinicalTrials.gov](https://ClinicalTrials.gov) under the reference code NCT02268695.

2.2 | Data

Data on 526 patients were gathered from patient records struc-
tured by a Case Report Form (CRF). Tumor response was
measured every 8 weeks after the start of the treatment. Each
measurement was performed by Computed Tomography (CT)
scan or MRI for neck and CT scan for chest, and abdomen until
disease progression. The sum of the largest diameters (SLDs) of
target lesions and the best overall response (BOR) were com-
puted according to RECIST v1.1. OS outcome was defined as
the time from the start of treatment to death or to the (right
censored) date of last follow-up and PPS as the duration from
the time of documented disease progression after the first line
treatment to death from any cause or to the (right censored)
date of last follow-up. Data about treatments (type of treatment,
treatment line, and toxicity events), disease characteristics
(histology, metastasis status, loco-regional recurrence, tumor
node metastasis (TNM) classification of the initial disease,
TNM classification of the loco-regional recurrence), clinical,
demographic, and epidemiological data (gender, performance
status, tobacco status, alcohol status, age, and country) were
also retrieved from the records. From the original set of over
40 baseline clinical and demographical features extracted from
the data, 12 were retained for analysis based on their clinical
relevance and their imbalance ratio, which refers to the extent
of uneven distribution among different categories within these
features. The dataset initially included seven patients who died
before receiving the first treatment (three from the training set
and four from the test set); these records were excluded from
the analysis. Following this, the dataset exhibited 58 missing
observations (1.32% missingness) in the training set and 23
missing observations (1.24% missingness) in the test set. To
ensure robust model training and evaluation, rows containing
missing values were also excluded. This step reduced the train-
ing set from 365 to 307 records and the test set from 154 to 131
records. Although the overall percentage of missingness was
relatively low, these exclusions were necessary to maintain the

integrity of the subsequent analyses and to avoid potential bi-
ases associated with imputation methods.

2.3 | Kinetics-Machine Learning

The predictive model was based on a specific two-step mod-
eling approach called the kinetics-Machine Learning [17–19]
(Figure 1B). It combines phenomenological modeling of the tumor
response and ML techniques for survival analysis to build pre-
dictive models. The time evolution of the SLDs (tumor kinetics,
TK) was modeled by tumor growth inhibition models [12, 13, 20].
Mixed-effects statistical learning techniques were used to esti-
mate the individual parameters. Finally, using different ML meth-
ods for survival analysis, the estimated tumor kinetic parameters
(patient-specific kinetic signature) were combined with relevant
baseline features to predict individual patients' survival.

2.4 | Tumor Kinetics Models

2.4.1 | Structural Models

Multiple TK models were assessed for their ability to fit the ob-
served SLD measurements. Let t, the time variable (months), TS,
the tumor size (i.e., SLD, mm), TS₀ (mm), the tumor size at the
first SLD assessment and τ (months), the time lag between first
SLD measurement and the treatment start (t = 0). In the first
model, tumor on-treatment regression or shrinkage is character-
ized by an exponential decrease with decay rate KS (months⁻¹)
and tumor growth is described by an exponential growth with
a growth rate KG (months⁻¹) (Table 1, model 1) [20]. The two
other variants of this model introduce one more parameter
[12, 13]. In the first one, the growth and shrinkage processes are
respectively modulated by the factors α and 1 − α, respectively
describing proportions of resistant and non-resistant tumor
cells (Table 1, model 2). In the third model, the decrease in the
treatment effect is modeled by a decay rate λ (Table 1, model
3) [12, 13]. For each model, an additional parameter, Time To
re-Growth (TTG), known for its ability to predict survival in
other tumor types, was also assessed and tested as a predictor
of OS [13].
# INTRODUCTION

Head and Neck Squamous Cell Carcinomas (HNSCCs) arise from the mucosal epithelium of the oral cavity, pharynx and the larynx, presenting a significant clinical challenge due to their aggressive nature and the associated morbidity. Regarding resectable HNSCCs, the standard of care is based on surgery followed by and radiotherapy +/- chemotherapy ¹,². However, for unresectable loco-regional HNSCCs, chemo-radiotherapy is the preferred modality. For patients with local or metastatic recurrence unable to benefit from local treatments such as surgery or radiotherapy, conventional treatment includes a triad of cisplatin and 5-fluorouracil (5-FU) combined with cetuximab (EXTREME regimen) ³ or a triad of cisplatin/carboplatin and 5-fluorouracil (5-FU) combined with pembrolizumab depending on the PD-L1 status ⁴. These regimens are typically dispensed over six cycles, each spanning 21 days, as a first line therapeutic approach. Moreover, a GORTEC 2014-01 TPExtreme clinical trial whose objective was to compare the efficacy and safety of the TPEx regimen with EXTREME regimen, could provide an alternative to standard of care in the first-line treatment of patients with recurrent or metastatic HNSCC, especially for those who might not be good candidates for up-front pembrolizumab treatment ⁵. Despite these therapeutic strides, the grim reality of a median survival oscillating between 12 to 15 months post-treatment⁵ underlines a pressing need for more precise prognostic models to guide clinical decisions and tailor therapeutic regimens.

Prediction of overall survival (OS) from early surrogate markers is of critical importance, both for clinical care and drug development (for instance to forecast the final issue of a phase 3 trial from early interim data). Current surrogate markers to predict OS are based on the response evaluation criteria in solid tumors (RECIST)⁶ that establish best overall response (BOR) using the relative change to baseline of the sum of the largest diameters (SLD) of target lesions. This also allows to define progression-free survival (PFS). However, such criteria have limitations. For instance, it has been observed that OS results could be uncorrelated to PFS or best overall response (BOR)⁷,⁸.

Furthermore, in the quest for enhanced prognostic prediction models, on-treatment longitudinal SLD measurements were used to model TK using semi-mechanistic “tumor growth inhibition” (TGI) models⁹,¹⁰. These models offer a glimpse into the dynamics of tumor response to treatment over time, thereby providing TK parameters for prognostic modeling. Multiple studies combining estimated TK parameters to survival analysis (TK–OS) revealed the importance of TK parameters in predicting OS, e.g., in metastatic colorectal cancers or non-small cell lung cancer (NSCLC)⁹,¹¹,¹². However, no direct comparison between the predictive abilities of TK parameters versus RECIST BOR has yet been performed.
<table><caption>**TABLE 1** | Different TK models.</caption><thead><tr><th>Name</th><th>Model</th><th>Parameters</th><th>TTG</th></tr></thead><tbody><tr><td>Model 1</td><td>TS(t, Θ) = TS<sub>0</sub>(e<sup>KG(t-τ)</sup> + e<sup>-KS(t-τ)</sup> - 1)</td><td>TS<sub>0</sub>, KG, KS</td><td>log(KS/KG) / (KS + KG) + τ</td></tr><tr><td>Model 2</td><td>TS(t, Θ) = TS<sub>0</sub>(αe<sup>KG(t-τ)</sup> + (1 - α)e<sup>-KS(t-τ)</sup>)</td><td>TS<sub>0</sub>, KG, KS, α</td><td>log(KS(1-α)/αKG) / (KS + KG) + τ</td></tr><tr><td>Model 3</td><td>TS(t, Θ) = TS<sub>0</sub>e<sup>(KS/λ)(1-e<sup>-λ(t-τ)</sup>)</sup> + KG(t-τ)</td><td>TS<sub>0</sub>, KG, KS, λ</td><td>log(KS/KG) / λ + τ</td></tr></tbody></table>

Abbreviation: TTG, Time To re-Growth.

## 2.4.2 | Statistical Model

The framework of the Nonlinear Mixed Effects Modeling (NLME) was used as a statistical tool to describe the inter-individual variability in the observed tumor kinetics. Let $Y^i = \{y_{i1}, \dots, y_{in}\}$, for $i \in \{1, \dots, N\}$, the vector of longitudinal SLD measurements for the $i$th patient, where $N$ represents the total number of patients and $n^i$, the number of measurements performed for patient $i$. The statistical model used in the fitting step was given by

$$y_{ij} = \mathrm{TS}^{(k)}(t_{ij}; \theta_{i}^{(k)}) + e_{ij}^{(k)} \quad (1)$$

where $\mathrm{TS}^{(k)}(t_{ij}, \theta_i)$ is the value of a specific structural model $k$ (for $k \in \{1, 2, 3\}$) (Table 1) at time $t_{ij}$, $\theta_i^{(k)} \in \mathbb{R}^p$ is the corresponding vector of parameters specific to the individual $i$ and $e_{ij}^{(k)}$ represents the residual error model. The vector of individual parameters $\theta_i^{(k)}$ was described by the combination of fixed effects $\theta_{pop}^{(k)}$, which are constant among the population and random effects $\eta_i^{(k)}$ which describe the inter-individual variability observed in the SLD kinetics. The random effects were assumed to be normally distributed with mean zero and variance-covariance matrix $\Omega^{(k)} = \mathrm{diag}(\omega_1, \dots, \omega_{l(k)})$ (where $l(k)$ refers to the number of kinetic parameters in a specific model $k$). Furthermore, the positivity of the parameters was ensured by assuming a log-normal distribution on the individual parameters. More precisely,

$$\log \theta_i^{(k)} = \log(\theta_{pop}^{(k)}) + \eta_i^{(k)}, \eta_i^{(k)} \sim \mathcal{N}(0, \Omega^{(k)}) \quad (2)$$

For the different structural models, the error model was assumed to be constant and was defined as follows:

$$e_{ij}^{(k)} = \sigma^{(k)} \varepsilon_{ij} \quad (3)$$

where $\varepsilon_{ij} \sim \mathcal{N}(0, 1)$ represents the residual error and $\sigma^{(k)}$ its standard deviation. In summary, the vector of population parameters is:

$$\psi^{(k)} = (\theta_{pop}^{(k)T}, \sigma^{(k)}, \omega^{(k)T})$$

where $\theta_{pop}^{(k)T}$ and $\omega^{(k)T}$ represent respectively the entries of the vector of fixed effects parameters and the entries of the vector of the random effects standard deviations.

## 2.4.3 | Model Calibration and Selection

The parameters $TS_0$ and $\tau$ were retrieved from the data. Some SLD measurements which were censored and set to zero due to the maximum CT (or MRI) scan slice thickness were reset to an offset of 2.5 mm (half of the value required for a lesion to be considered as measurable according to RECIST v1.1) and marked as left-censored. The contribution of censored measurements to the observed likelihood function was handled by lower Limit Of Quantification (LOQ) censoring [21]. Furthermore, the population parameters $\psi^{(k)}$ of each statistical model were estimated by computing the maximum likelihood estimate $\hat{\psi}^{(k)}$ of the observed likelihood function with the Stochastic Approximation of the Expectation-Maximization Algorithm implemented in the Monolix R API (Version 2021R2). The Fisher information matrix was derived from the observed likelihood by a Markov Chain Monte Carlo algorithm implemented in Monolix. The relative standard errors (RSE) and corresponding confidence intervals were computed for each component of $\hat{\psi}^{(k)}$. The estimator $\hat{\theta}_i^{(k)}$ of the individual kinetic parameters was defined as the mode of the posterior conditional distribution $p(\theta_i^k | y_{ij}, \hat{\psi}^{(k)})$.

Goodness-of-fit of the NLME model was assessed by the Bayesian information criterion (BIC). A RSE of $\le 50\%$ was established as a threshold to ensure parameter stability and model identifiability, thereby indicating that the parameter estimates were reliable and precise. Additionally, diagnostic plots, including individual curve fits and comparisons of observed versus predicted SLD measurements, were employed to illustrate and validate the model's performance.

## 2.5 | ML for Survival Analysis

The data on the 526 patients were split into a training and test set, which represented respectively 70% and 30% of the full dataset. Initial data pre-processing was performed on the training set, from which the scaling parameters—specifically the mean and standard deviation of continuous variables—were computed to standardize the features. These scaling parameters were then applied to the test set to maintain consistency. TK model selection and multiple survival ML algorithms were benchmarked on the training set. The best TK model selected on the training set was calibrated on the test set using the NLME priors obtained on the training set and the corresponding kinetic parameters were retrieved to assess the accuracy of the best ML model selected
In addition, a slew of biases threatens the improvements afforded by the TK parameters on predicting OS. One such bias is the time-dependent covariate bias <sup>13,14</sup>, which can skew the derived insights <sup>15</sup>. Amidst these backdrops, post-progression survival (PPS) arises as an adequate prognostic outcome. The pivotal rationale favoring PPS over OS is mainly anchored in the attempt to mitigate the confounding impact that subsequent lines of therapy can have on OS <sup>7</sup>. Unlike PFS, PPS is strongly correlated to OS (r = 0.97, p < 0.05, R<sup>2</sup> = 0.94) <sup>8</sup>. For instance, in cases of NSCLC, an increment in PFS does not unequivocally translate to an enhancement in OS, however, PPS exhibits a strong association with OS following the first-line treatment <sup>8</sup>. Moreover, PPS has been acknowledged for its potential in alleviating the immortal time bias associated with the employment of model-derived metrics as covariates for survival analysis <sup>15</sup>.

This is mostly because such approaches where initially developed for drug development, in which case such bias does not apply. This is mainly because the output of interest was at the study or treatment arm level, rather than the individual one. For individual survival predictions, one needs to account for the time-dependent covariate bias that arises when *on-treatment* longitudinal data (i.e., TK) is used to predict at a *pre-treatment* time point, which is what is done when the outcome is OS <sup>15</sup>. Such approach is mathematically ill-defined as it uses data from the future to predict in the present. In the context of TK-OS, using all the available longitudinal data for predicting OS is flawed because already the total duration of the "on-treatment" phase is a strong predictor of OS. To circumvent this, our approach entails the estimation of TK parameters from SLD measurements restricted to a specific duration and thereafter predicting the subsequent survival outcome. Specifically, we examined two well-posed operational prediction problems: 1) forecasting post-cycle 4 OS (OS4) from tumor kinetics data up to that point (TK4 – OS4) and 2) predicting post-progression survival (PPS) from first line tumor kinetics (TK1L – PPS).

To address these, we harnessed the strengths of population TK modeling, utilizing non-linear mixed effects (NLME), coupled with machine learning (ML) survival analysis. After benchmarking the currently available ML survival models, we compared the predictive power of TK parameters versus Best Overall Response (BOR) for survival prediction.
on the training set. The estimated individual kinetic parame-
ters—including TTG—were added to the baseline features for
the assessment of different predictive ML and statistical models
for survival analysis.

To avoid skewness and to be consistent with the log-normal hy-
pothesis in the modeling process, a logarithmic transformation
was applied to the parameters KG and KS.

**2.5.1 Benchmarking Procedure**

Eight predictive ML models—CoxTime [22], DeepHit [23], DeepSurv [24], LogHaz [25, 26], PC-Hazard [27], Survival SVM [28, 29], XGBoost [30], and Random Survival Forest (RSF) [31]—and five classical parametric and semi-parametric survival models—Cox Proportional Hazard (PH), Accelerated Failure Time (AFT) [32], Cox Lasso [33], CoxBoost, and GLMBoost [34] and [35]—were assessed for their ability to predict PPS. A brief description of each model and its methodological approach is provided in Table S1. The aforementioned models were compared by computing Harrell's Concordance index (c-index) [36] using nested and repeated k-fold cross-validation. Each inner loop contained 30 repeats of twofold cross-validations, used to tune the ML models. Each outer loop consisted of 10 repeats of fivefold cross-validations, used to assess the generalizability of the different models. The models were tuned by performing a random search on the numerical hyperparameters domain. Overall, 50 training and test sets were used to compute the predictive performance of each model. The benchmarking procedure and the different models were implemented using the *mlr3* packages in *R* version 4.1.0.

To assess the features importance, a permutation-based method was applied to the best ML model on the training set. However, for cases where Cox PH emerged as the best model, the Hazard Ratios (HR) were used to quantify feature significance. To stratify the survival curves based on TK parameters, the maximally selected rank statistic from the *maxstat R* package was employed.

**2.5.2** | Test Step

First, the best TGI model obtained on the training set was cal-
ibrated on the test set using the priors obtained on the train-
ing set and the estimated kinetic parameters were added to the
baseline features. Afterwards, using 1000 bootstrap sets derived
from the test set, the accuracy of the best performing ML model
was assessed on the test set by computing the c-index and its 95%
confidence interval.

**2.5.3 | Comparison of Different Predictor Models**

The best model selected from the benchmarking step was used
to assess the predictive quality of different sets of predictors:
baseline (BSL) features, BSL+BOR, BSL+TK, and TK. These
models were assessed on the training set using 10 repeats of
fivefold cross validations. The generalization quality of the
models was assessed on the test set using 1000 bootstrap sets.

Furthermore, a Kruskal–Wallis test was used to compare the c-
index of the models.

**2.6 | Unbiased Two-Step Modeling Approaches: TK4-OS4 and TKL1-PPS**

It would be biased to predict the time-to-death of a patient from baseline (i.e., OS) using TK data that is not available at this time. To avoid this, we introduced two unbiased two-step modeling approaches. The first, denoted TK4-OS4, used 4 cycles SLD data to estimate TK parameters (TK4) and predicts post cycle 4 OS (OS4). The second approach, denoted TKL1-PPS, used TK parameters identified from longitudinal SLD measurements restricted to the first line of treatment + the maintenance period (e.g., the time to disease progression), to predict PPS.

3 | Results

**3.1 | Patients and Disease Characteristics**

Overall, 526 patients treated with either the EXTREME or TPEx regimen were retrieved from the GORTEC 2014-01 TPExtreme clinical trial [5]. Patients and disease characteristics are summarized in Table S3. In this study, 60.1% of the patients received a second-line treatment while 33.2% received a third-line treatment, and finally 13.3% received a fourth-line treatment. In this study, 44.6% of the primary tumors were in the oropharynx and 64% of the patients had metastasis at inclusion (Table S3). The typical treatment history of a patient is depicted in Figure 1A. It highlights the inter-patient heterogeneity in terms of both efficacy and toxicity, thus emphasizing the complexity of clinical management of such patients at bedside.

**3.2 | Tumor Kinetics Modeling**

The SLD measurements on 365 (70%) patients—restricted to four cycles (total number of time points = 672 and median number of time points per patient = 2)—were used to calibrate the different TK models. Among the tested models, Model 1 demonstrated the best agreement with the data, as indicated by the lowest values for all information criteria (AIC, BIC, and BICc), suggesting a superior balance between model fit and complexity. Furthermore, model 1 was found to be the best to describe the data in the different cases (TK4 and TKL1) (Table S1 for TK4). All the parameters were identifiable with a maximum RSE of 8.73%, confirming model stability and parameter identifiability (Table S2).

Diagnostic plots did not reveal any misspecification of the
model at the individual level, supporting the adequacy of model
1 for describing the observed tumor dynamics (Figure 2A,B).
Recalibrating the model using the treatment arm as a covariate
to the kinetic parameters (KS, KG) revealed a potential effect of
the TPEx treatment on the tumor growth rate (KG) in compari-
son to its effect on the tumor shrinkage rate (KS). Indeed, a Wald
test demonstrated a statistically significant association between
the treatment arm and KG, (p = 6.35 × 10⁻³) But no association
was found between the treatment arm and KS. Furthermore,
# MATERIALS AND METHODS

## Patients

Data consisted of all patients enrolled in the GORTEC 2014-01 TPExtreme clinical trial ⁵. Included patients had at least one measurable lesion according to the Response Evaluation Criteria in Solid Tumors (RECIST) v1.1 ⁶, ECOG performance status of 1 or 0, and were aged 18 – 70 years with histologically confirmed metastatic or recurrent squamous-cell carcinoma of the oral cavity, hypopharynx, oropharynx, or larynx, unsuitable for loco-regional curative treatments. A total of 526 participants were randomly assigned (1:1) to each treatment group. On one hand, the TPEx regimen was administrated and consisted of 4 cycles of a combination of docetaxel, cisplatin, and cetuximab. On the other hand, the EXTREME regimen consisted of 6 cycles of the standard of care (fluorouracil, cisplatin and cetuximab) (Figure 1.A).

All enrolled patients gave written informed consent before any study procedure. The trial was carried out in respect of good clinical practice guidelines and the declaration of Helsinki. It was approved by competent authorities and ethics committees in France, Spain, and Germany and registered with ClinicalTrials.gov under the reference code NCT02268695.

## Data

Data on 526 patients were gathered from raw patients' records structured by a Case Report Form (CRF). Tumor response was measured every 8 weeks after the start of the treatment. Each measurement was performed by Computed Tomography (CT) scan or MRI for neck and CT scan for chest, and abdomen until disease progression. The sum of largest diameters (SLDs) of target lesions and best overall response (BOR) were computed according to RECIST v1.1. OS outcome was defined as the time from the start of treatment to death or to the (right censored) date of last follow-up and PPS as the duration from the time of documented disease progression after the first line treatment to death from any cause or to the (right censored) date of last follow-up. Data about treatments (type of treatment, treatment line, toxicity events), disease characteristics (histology, metastasis status, loco-regional recurrence, tumor node metastasis (TNM) classification of the initial disease, TNM classification of the loco-regional recurrence), clinical, demographic, and epidemiological data (gender, performance status, tobacco status, alcohol status, age, country) were also retrieved from the raw records. Over 40 baseline clinical and demographical features gathered from the raw data, 15 were retained for analysis considering their clinical relevance and their imbalance ratio (IR) (Error! Reference source not found.). The total fraction of missingness in the dataset was 1.24% (1.45% on the test set). Therefore, due to their small amount, missing values were discarded from the training and the test set.
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Diagnostic plots (A) Representative individual fits. The purple line represents the model prediction and the blue dots the observed SLDs. The red, green and dark-blue dotted vertical lines represent the end of the first, second line and third lines, respectively. (B) Observed SLDs vs. predicted SLDs.</figcaption></figure>

incorporating the treatment arm or any additional covariates to the model resulted in a decreased goodness-of-fit, introduced iden- tifiability issues and did not show any significant improvement in the survival analysis. Given the goal of predicting outcomes on the entire patient cohort without introducing covariate-related bias, a model without covariates was kept for the predictive analysis. This choice ensured a more generalizable model for the popula- tion, regardless of specific treatment arm effects. The differences in the models did not significantly alter the overall conclusions regarding TK parameters and survival predictions.

Figure 3 and 0.601, 0.598, and 0.595 for TKL1-OS, Figure S1). However, an examination of the results on the test set for TK4-OS4 revealed that RSF exhibited the best generalization properties in comparison to Cox PH and AFT (Figure S2). Overall, TKL1-PPS (best training c-index=0.6) exhibited worse performances than TK4-OS4 (best training c-index=0.65, Figure 3B). Moreover, deep learning architectures, specifically DeepHit and DeepSurv, exhib- ited suboptimal performances for both TK4-OS4 and TKL1-PPS. Upon examination of the test set, for TK4-OS4, the Cox PH model achieved a mean c-index of 0.54 (95% CI: 0.52–0.56) (Figure S2).

## 3.3 | ML For Survival Analysis

### 3.3.1 | ML Models Benchmarking Unbiased Predictions Using TK4-OS4 and TKL1-PPS

To avoid time-dependent covariate bias, we considered either predicting post-cycle 4 survival (OS4) from data restricted to the first four treatment cycles (TK4-OS4) or predict post-progression survival (PPS) from TK parameters derived from the data during the first line of treatment until progression (TKL1-PPS). On the training set, for both TK4-OS4 and TKL1-OS, the random sur- vival forest (RSF) model did not significantly outperformed Cox PH or AFT (median c-indices = 0.65, 0.635, and 0.63 for TK4-OS4,

### 3.3.2 | Feature Importance

Noteworthy, the most important feature for OS4 prediction was first the patient performance status followed by the TK param- eters (logKG, logKS, and TTG) and the baseline tumor size TS₀ (Figure 4A). In addition, utilizing the maximally selected rank statistic, TS₀, logKG and logKS provided a significant separation of the survival curves (log-rank *p* < 0.01 and *p* < 0.05, respec- tively) (Figure 4B–D).

However, for TKL1-PPS, the top 3 parameters were loco-regional relapse, followed by TS₀, the performance status (Figure S3A). Surprisingly, Unlike TK4-OS4, logKS, triggered crossed survival
## Kinetics-Machine Learning

The predictive model was based on a specific two-step modeling approach called the kinetics-Machine Learning <sup>16-18</sup> (Figure 1.B). It combines phenomenological modeling of the tumor response and ML techniques for survival analysis to build predictive models. The time evolution of the SLDs (tumor kinetics, TK) was modeled by tumor growth inhibition models <sup>19,9,10</sup>. Mixed-effects statistical learning techniques were used to estimate the individual parameters. Finally, using different ML methods for survival analysis, the estimated tumor kinetic parameters (patient-specific kinetic signature) were combined with relevant baseline features to predict individual patients' survival.

A

<figure><img src="image_2.png" /></figure>

B

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Treatment history and kinetics-ML prediction of survival. (A) Treatment history of a sample of patients. Cetuximab was the treatment given during maintenance. (B) Steps describing mechanistic learning. It combines modeling tumor kinetics and ML for survival prediction.</figcaption></figure>

## Tumor kinetics models

### Structural models

Multiple TK models were assessed for their ability to fit the observed SLD measurements. Let $t$, the time variable (*months*), $TS$, the tumor size (i.e., SLD, *mm*), $TS_0$ (*mm*), the tumor size at the first SLD assessment and $\tau$ (*months*), the time lag between first SLD measurement and the treatment start ($t = 0$). In the first model, tumor on-treatment regression or shrinkage is characterized by an exponential decrease with decay rate $KS$ (*months*⁻¹) and tumor growth is described by an exponential growth with a growth rate $KG$ (*months*⁻¹) (*Table 1*, model 1) <sup>19</sup>. The two other variants of this model introduce one more parameter <sup>9,10</sup>. In the first one, the growth and shrinkage processes are respectively modulated by the factors $\alpha$ and $1 - \alpha$,
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Benchmark of ML methods, TK4-OS4. (A) Graphical representation of the different ML and statistical methods. The green boxes are the tested models. (B) Boxplots of the results of the benchmarking process on OS4 prediction. The red Diamond is the median c-index and the middle line the mean c-index.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

curves (Figure S3D) potentially indicating a time-dependent na-
ture of the growth rate parameter.

**3.3.3 | Comparison of Different Predictor Models: TK4 Parameters Are Better Predictors of OS4 Than BOR**

To compare the relative predictive performances of BOR and TK4 as predictors of OS4, four sets of predictors were considered in combination with RSF: 12 baseline features (BSL, see Table S4) alone, BSL+BOR, BSL+TK4 and TK4 (Figure 5). In the training set, the BSL+TK4 outperformed BOR+BSL 0.63 (95%CI: 0.61 – 0.65) versus 0.612 (95%CI:0.58 – 0.64), respectively, *p* < 0.001, Kruskal-Wallis, (Figure 5A). This trend persisted in the test dataset, with respective c-index of 0.63 (95% CI: 0.61 – 0.65) and 0.52 (95% CI: 0.50 – 0.54) (*p* < 0.0001, Figure 5B). The integration of BSL with BOR adversely impacted the model accuracy on test set as compared to BSL alone (Figure 5B). However, this was not observed in TKL1-PPS, such behavior was absent (Figure S4A,B).

4 | Discussion

Tumor kinetics modeling to predict overall survival (TK-OS) has
been extensively developed in the last 15 years [37]. Multiple can-
cer types, often in the advanced or metastatic stage have been
studied, including prostate [38], breast [39], colorectal [12], non-
small cell lung cancer (NSCLC) [40], or renal cell carcinoma [41].
Our study is the first to apply such modeling to HNSCC, using

data from the TPExtreme trial [5]. In addition, it is one of the few
to consider and benchmark machine learning survival models
versus classical (semi-)parametric survival analysis methods [42].

As often observed for the other cancer types, such as NSCLC [14], breast, and colorectal cancer [39], on-treatment TK could be accurately described by a simple pattern mathematically: the sum of two exponentials. A decreasing phase is governed by a shrinkage parameter KS and a regrowth phase is led by a parameter KG.

The primary objective of the TK modeling was to assess the predictive value of the TK parameter on the survival outcomes. Incorporating the treatment arm as a covariate to investigate its potential impact on tumor dynamics resulted in a reduced model performance, characterized by decreased goodness-of-fit and identifiability issues. Consequently, covariates (including treatment arms) were excluded from the final model to ensure identifiability and robustness. Despite these limitations, it is noteworthy that the exploratory analysis of treatment arms as covariates using simulated SLDs revealed a significant effect of the TPEx treatment on the tumor growth rate (KG), indicating a potential differential impact of the treatment regimens on tumor dynamics. Specifically, Using the treatment arm as a covariate demonstrated a significantly slower regrowth (parameter KG) and deeper SLD decrease in TPEx arm versus the control EXTREME, particularly within the first 3 months. Additionally, observations of simulated median TK suggested a benefit of TPEx versus EXTREME, observed at ~12 months, with a better safety profile [5].

Unlike several previous TK modeling studies mostly aimed
to assist drug development [12, 14, 39, 40, 43], here we turned
respectively describing proportions of resistant and non-resistant tumor cells (Table 1, model 2). In the third model, the decrease in the treatment effect is modeled by a decay rate $\lambda$ (Table 1, model 3) ⁹,¹⁰. For each model, an additional parameter, Time To re-Growth (TTG), known for its ability to predict survival in other tumor types, was also assessed and tested as a predictor of OS ¹⁰.

<table><caption>Table 1. Different TK models</caption><thead><tr><th>Name</th><th>Model</th><th>Parameters</th><th>TTG<sup>1</sup></th></tr></thead><tbody><tr><td>Model 1</td><td>TS(t, Θ) = TS<sub>0</sub>(e<sup>KG(t-τ)</sup> + e<sup>-KS(t-τ)</sup> - 1)</td><td>TS<sub>0</sub>, KG, KS</td><td>log(<sup>KS</sup>/<sub>KG</sub>) / (KS + KG) + τ</td></tr><tr><td>Model 2</td><td>TS(t, Θ) = TS<sub>0</sub>(αe<sup>KG(t-τ)</sup> + (1 - α)e<sup>-KS(t-τ)</sup>)</td><td>TS<sub>0</sub>, KG, KS, α</td><td>log(<sup>KS(1 - α)</sup>/<sub>αKG</sub>) / (KS + KG) + τ</td></tr><tr><td>Model 3</td><td>TS(t, Θ) = TS<sub>0</sub>e<sup>-(KS/λ)(1 - e<sup>-λ(t-τ)</sup>) + KG(t-τ)</sup></td><td>TS<sub>0</sub>, KG, KS, λ</td><td>log(<sup>KS</sup>/<sub>KG</sub>) / λ + τ</td></tr></tbody></table>

## Statistical model

The framework of the Nonlinear Mixed Effects Modelling (NLME) was used as a statistical tool to describe the inter-individual variability in the observed tumor kinetics. Let $Y^i = \{y_{i1}, \dots, y_{in^i}\}$, for $i \in \{1, \dots, N\}$, the vector of longitudinal SLD measurements for the $i$th patient, where $N$ represents the total number of patients and $n^i$, the number of measurements performed for patient $i$. The statistical model used in the fitting step was given by

$$y_{ij} = TS^{(k)}(t_{ij}; \theta_i^{(k)}) + e_{ij}^{(k)}, \quad (1)$$

where $TS^{(k)}(t_{ij}, \theta_i)$ is the value of a specific structural model $k$ (for $k \in \{1, 2, 3\}$)(Table 1) at time $t_{ij}$, $\theta_i^{(k)} \in \mathbb{R}^p$ is the corresponding vector of parameters specific to the individual $i$ and $e_{ij}^{(k)}$ represents the residual error model. The vector of individual parameters $\theta_i^{(k)}$ was described by the combination of fixed effects $\theta_{pop}^{(k)}$, which are constant among the population and random effects $\eta_i^{(k)}$ which describe the inter-individual variability observed in the SLD kinetics. The random effects were assumed to be normally distributed with mean zero and variance-covariance matrix $\Omega^{(k)} = \text{diag}(\omega_1, \dots, \omega_{l^{(k)}})$ (where $l^{(k)}$ refers to the number of kinetic

<footnote n="1">¹ Time To re-Growth</footnote>
parameters in a specific model k). Furthermore, the positivity of the parameters was ensured by assuming a log-normal distribution on the individual parameters. More precisely,

$$ \log \theta_i^{(k)} = \log(\theta_{pop}^{(k)}) + \eta_i^{(k)}, \quad \eta_i^{(k)} \sim \mathcal{N}(0, \Omega^{(k)}). \quad (2) $$

For the different structural models, the error model was assumed to be constant and was defined as follows:

$$ e_{ij}^{(k)} = \sigma^{(k)} \varepsilon_{ij}, \quad (3) $$

where $\varepsilon_{ij} \sim \mathcal{N}(0, 1)$ represents the residual error and $\sigma^{(k)}$ its standard deviation. In summary, the vector of population parameters is:

$$ \psi^{(k)} = (\theta_{pop}^{(k)T}, \sigma^{(k)}, \omega^{(k)T}) $$

where $\theta_{pop}^{(k)T}$ and $\omega^{(k)T}$ represent respectively the entries of the vector of fixed effects parameters and the entries of the vector of the random effects standard deviations.

## Model calibration and selection

The parameters $TS_0$ and $\tau$ were retrieved from the data. Some SLD measurements which were censored and set to zero due to the maximum CT (or MRI) scan slice thickness were reset to an offset of 2.5 mm (half of the value required for a lesion to be considered as measurable according to RECIST v1.1) and marked as left-censored. The contribution of censored measurements to the observed likelihood function was handled by lower Limit Of Quantification (LOQ) censoring $^{20}$. Furthermore, the population parameters $\psi^{(k)}$ of each statistical model were estimated by computing the maximum likelihood estimate $\hat{\psi}^{(k)}$ of the observed likelihood function with the Stochastic Approximation of the Expectation-Maximization Algorithm implemented in the Monolix R API (Version 2021R2). The Fisher information matrix was derived from the observed likelihood by a Markov Chain Monte Carlo algorithm implemented in Monolix. The relative standard errors and corresponding confidence intervals were computed for each component of $\hat{\psi}^{(k)}$. The estimator $\hat{\theta}_i^{(k)}$ of the individual kinetic parameters was defined as the mode of the posterior conditional distribution $p(\theta_i^k | y_{ij}; \hat{\psi}^{(k)})$.

Goodness-of-fit of the NLME model was assessed by the Bayesian information criterion and the parameters were considered identifiable when their RSE did not exceed 50 %.

## ML for survival analysis

The data on the 526 patients were split into a training and test set, which represented respectively 70% and 30% of the full dataset. Initial data pre-processing was performed on the training set, from which the scaling parameters were extracted. These parameters were then
<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | Comparison of different predictive models for TK4-OS4. The red diamond and middle line are the median and mean c-index, respectively. (A) Predictive models on the training set. (B) Predictive models on the test set. *p* refers to the *p*-value of the Kruskal-Wallis test.</figcaption></figure>

survival chances of a patient are predicted to be very small, the clinician or the patient might decide to stop anti-cancer therapy. However, to really decipher the impact of the first line TK on the response to second-line options would require the knowledge of the treatments received post-first line. A limitation of our study is that such data was not available in the trial records. Further investigations should be conducted to address this question.

To maximize predictive performances and account for possible nonlinear interactions between predictors of survival, we benchmarked the current state-of-the art survival machine learning for the two survival prediction problems. We compared traditional (Cox proportional hazard) modeling to tree ensemble models (e.g., gradient boosting or random survival forest) and neural networks (CoxTime, DeepSurv, and DeepHit). Our results indicate that the more complex (and less interpretable) models only rarely and marginally surpassed the standard, linear (Cox) model. This contrasts with another work with more patients and variables, where a random survival forest outperformed a Cox model [17]. However, they align with other studies showing that complex ML might not overcome classical statistical models for cases relative to predictive modeling in healthcare from structured data, especially when the number of samples is small (of the order of 100–1000) [47–49]. Rather, advanced neural network engineering performs better for computer vision or natural language processing tasks [49]. Yet, the best model for TK4-OS4 was a random survival forest model. It demonstrated a modest (but significant) predictive performance of individual survival with a c-index of 0.592. It is important to note that the relatively modest c-index reflects the inherent complexity of HNSCC and the heterogeneity of the patient population. These findings underscore the necessity of applying the mechanistic learning approach to a larger and more diverse dataset, incorporating additional biomarkers.

When reducing the number of variables to only BSL + TK4, the model: (1) did not exhibit signs of overfitting (same performances in the training and test sets) but most importantly (2) indicated a substantially better prediction of OS4 than a model using BSL + BOR. This last point indicates that the current standard RECIST guidelines for evaluation and monitoring of response could be improved using dynamic TK-based metrics that account for the entire on-treatment longitudinal data.

Unfortunately, our results indicated a limited value of TKL1 model-based metrics for PPS prediction. However, several refinements could be further investigated. For example, modeling individual lesion sizes rather than the SLD unraveled substantial inter-lesion and inter-organ variability within a specific patient, and this variability was linked to response [50–52]. In addition, considering the total volume rather than one-dimensional measurements can reveal inter-lesion variability that could not be observed with diameters [53] and provide a better assessment of the treatment effect [54]. Spatiotemporal heterogeneity across metastases has also been found relevant in a large study comprising more than 11,000 lesions in 2802 colorectal cancer patients [55]. Further, modeling the kinetics of functional tumor markers could add predictive value. For instance, Schindler et al. used the per-lesion maximal standard uptake value (SUV<sub>max</sub>) obtained from positron emission tomography and not only found significant inter-lesion variability but also that this marker outperformed tumor size for prediction of survival [50].

Beyond invasive or irradiating tumor markers, circulating markers could also enhance predictions [56–59]. We recently found improved survival predictive performances when adding four simple blood markers kinetics (BK) from routine lab tests to TK data in a large dataset of non-small cell lung cancer patients treated with an immune-checkpoint inhibitor [17, 19]. The data comprised three phase 2 trials used as the training set for model development (862 patients) and a phase 3 trial used as the test set for model evaluation (377 and 354 patients in the investigational and control (docetaxel) arms, respectively). One of the main advantages of such BKs is that they are non-invasive and can be sampled more frequently than tumor measurements. In addition, the advent of liquid biopsies provides qualitative and quantitative data on circulating free or tumor DNA that pave a promising avenue to develop more accurate predictive computational models leveraging on-treatment data [60, 61].

In conclusion, our work suggests that early analysis of first-line tumor kinetics through mechanistic learning provides valuable prognostic insights. Future studies should aim to integrate additional biomarkers and explore the utility of larger datasets to improve prediction accuracy and better guide treatment decisions.
utilized to standardize the test set. TK model selection and multiple survival ML algorithms
were benchmarked on the training set. The best TK model selected on the training set was
calibrated on the test set using the NLME priors obtained on the training set and the
corresponding kinetic parameters were retrieved to assess the accuracy of the best ML model
selected on the training set. The estimated individual kinetic parameters – including *TTG* –
were added to the baseline features for the assessment of different predictive ML and
statistical models for survival analysis.

To avoid skewness and to be consistent with the log-normal hypothesis in the modeling
process, a logarithmic transformation was applied to the parameters *KG* and *KS*.

## Benchmarking procedure

Eight predictive ML models – CoxTime ²¹, DeepHit ²², DeepSurv ²³, LogHaz ²⁴,²⁵, PC-Hazard ²⁵, Survival SVM ²⁶,²⁷, XGBoost and Random Survival Forest (RSF) ²⁸ – and five classical parametric and semi-parametric survival models – Cox Proportional Hazard (PH), Accelerated Failure Time (AFT), Cox Lasso, CoxBoost and GLMBoost ²⁹ and ³⁰ – were assessed for their ability to predict PPS. The aforementioned models were compared by computing Harrell’s Concordance index (c-index) ³¹ using nested and repeated k-fold cross-validation. Each inner loop contained 30 repeats of 2-fold cross-validations, used to tune the ML models. Each outer loop consisted of 10 repeats of 5-fold cross validations, used to assess the generalizability of the different models. The models were tuned by performing a random search on the numerical hyperparameters domain. Overall, 50 training and test sets were used to compute the predictive performance of each model. The benchmarking procedure and the different models were implemented using the *mlr3* packages in *R* version 4.1.0.

To assess the features importance, a permutation-based method was applied to the best ML
model on the training set. However, for cases were Cox PH emerged as the best model, the
Hazard Ratios (HR) were used to quantify feature significance. To stratify the survival curves
based on TK parameters, the maximally selected rank statistic from the *maxstat R* package
was employed.

## Test step

First, the best TGI model obtained on the training set was calibrated on the test set using the priors obtained on the training set and the estimated kinetic parameters were added to the baseline features. Afterwards, using 100 bootstrap sets derived from the test set, the accuracy of the best performing ML model was assessed on the test set by computing the c-index and its 95% confidence interval.
## Author Contributions

K.A., A.A., J.G., S.S., and S.B. wrote the manuscript. K.A., S.S., and S.B. designed the research. K.A. performed the research. K.A., S.B., S.S., and A.A. analyzed the data. K.A. contributed new reagents/analytical tools.

## Acknowledgments

This work is part of the QUANTIC project funded by ITMO Cancer AVIESAN and the French Institute National du Cancer (grant #19CM148-00). This work is part of the DIGPHAT project which was supported by a grant by a grant from the French government, managed by the National Research Agency (ANR), under the France 2030 program, reference ANR-22-PESN-0017.

## Conflicts of Interest

The authors declare no conflicts of interest.

## References

1. D. E. Johnson, B. Burtness, C. R. Leemans, V. W. Y. Lui, J. E. Bauman, and J. R. Grandis, "Head and Neck Squamous Cell Carcinoma," *Nature Reviews. Disease Primers* 6 (2020): 92.

2. S. Okano, A. Homma, N. Kiyota, et al., "Induction Chemotherapy in Locally Advanced Squamous Cell Carcinoma of the Head and Neck," *Japanese Journal of Clinical Oncology* 51 (2021): 173-179.

3. J. B. Vermorken, R. Mesia, F. Rivera, et al., "Platinum-Based Chemotherapy Plus Cetuximab in Head and Neck Cancer," *New England Journal of Medicine* 359 (2008): 1116-1127.

4. B. Burtness, K. J. Harrington, R. Greil, et al., "Pembrolizumab Alone or With Chemotherapy Versus Cetuximab With Chemotherapy for Recurrent or Metastatic Squamous Cell Carcinoma of the Head and Neck (KEYNOTE-048): A Randomised, Open-Label, Phase 3 Study," *Lancet* 394 (2019): 1915-1928.

5. J. Guigay, A. Aupérin, J. Fayette, et al., "Cetuximab, Docetaxel, and Cisplatin Versus Platinum, Fluorouracil, and Cetuximab as First-Line Treatment in Patients With Recurrent or Metastatic Head and Neck Squamous-Cell Carcinoma (GORTEC 2014-01 TPExtreme): A Multicentre, Open-Label, Randomised, Phase 2 Trial," *Lancet Oncology* 22 (2021): 463-475.

6. E. A. Eisenhauer, P. Therasse, J. Bogaerts, et al., "New Response Evaluation Criteria in Solid Tumours: Revised RECIST Guideline (Version 1.1)," *European Journal of Cancer* 45 (2009): 228-247.

7. Z. Zhang, C. Xie, T. Gao, Y. Yang, Y. Yang, and L. Zhao, "Identification on Surrogating Overall Survival With Progression-Free Survival of First-Line Immunochemotherapy in Advanced Esophageal Squamous Cell Carcinoma—An Exploration of Surrogate Endpoint," *BMC Cancer* 23 (2023): 145.

8. A. Shahnam, N. Hitchen, U. Nindra, et al., "Objective Response Rate and Progression-Free Survival as Surrogates for Overall Survival Treatment Effect: A Meta-Analysis Across Diverse Tumour Groups and Contemporary Therapies," *European Journal of Cancer* 198 (2024): 113503.

9. L. Belin, A. Tan, Y. De Rycke, and A. Dechartres, "Progression-Free Survival as a Surrogate for Overall Survival in Oncology Trials: A Methodological Systematic Review," *British Journal of Cancer* 122 (2020): 1707-1714.

10. E. Saad, A. Katz, K. Machado, and M. Buyse, "Post-Progression Survival (PPS) and Overall Survival (OS) According to Treatment Type in Contemporary Phase III Trials in Advanced Breast Cancer (ABC)," *Cancer Research* 69 (2009): 5116.

11. H. Imai, K. Mori, K. Wakuda, et al., "Progression-Free Survival, Post-Progression Survival, and Tumor Response as Surrogate Markers

for Overall Survival in Patients With Extensive Small Cell Lung Cancer," *Annals of Thoracic Medicine* 10 (2015): 61-66.

12. L. Claret, P. Girard, P. M. Hoff, et al., "Model-Based Prediction of Phase III Overall Survival in Colorectal Cancer on the Basis of Phase II Tumor Dynamics," *Journal of Clinical Oncology* 27 (2009): 4103-4108.

13. L. Claret, M. Gupta, K. Han, et al., "Evaluation of Tumor-Size Response Metrics to Predict Overall Survival in Western and Chinese Patients With First-Line Metastatic Colorectal Cancer," *Journal of Clinical Oncology* 31 (2013): 2110-2114.

14. L. Claret, J. Y. Jin, C. Ferté, et al., "A Model of Overall Survival Predicts Treatment Outcomes With Atezolizumab Versus Chemotherapy in Non-Small Cell Lung Cancer Based on Early Tumor Kinetics," *Clinical Cancer Research* 24 (2018): 3292-3298.

15. C. van Walraven, D. Davis, A. J. Forster, and G. A. Wells, "Time-Dependent Bias Was Common in Survival Analyses Published in Leading Clinical Journals," *Journal of Clinical Epidemiology* 57 (2004): 672-682.

16. H. Mistry, "Time-Dependent Bias of Tumor Growth Rate and Time to Tumor Regrowth," *CPT: Pharmacometrics & Systems Pharmacology* 5 (2016): 587.

17. S. Benzekry, M. Karlsen, C. Bigarré, et al., "Predicting Survival in Patients With Advanced NSCLC Treated With Atezolizumab Using Pre- and On-Treatment Prognostic Biomarkers," *Clinical Pharmacology & Therapeutics* 116 (2024): 1110-1120.

18. S. Benzekry, "Artificial Intelligence and Mechanistic Modeling for Clinical Decision Making in Oncology," *Clinical Pharmacology and Therapeutics* 108 (2020): 471-486.

19. S. Benzekry, M. Karlsen, A. El Kaoutari, et al., "Supporting Decision Making and Early Prediction of Survival for Oncology Drug Development Using a Pharmacometrics-Machine Learning Based Model," *PAGE* 30 (2022): 10276.

20. W. D. Stein, J. Yang, S. E. Bates, and T. Fojo, "Bevacizumab Reduces the Growth Rate Constants of Renal Carcinomas: A Novel Algorithm Suggests Early Discontinuation of Bevacizumab Resulted in a Lack of Survival Advantage," *Oncologist* 13 (2008): 1055-1062.

21. A. Samson, M. Lavielle, and F. Mentré, "Extension of the SAEM Algorithm to Left-Censored Data in Nonlinear Mixed-Effects Model: Application to HIV Dynamics Model," *Computational Statistics & Data Analysis* 51 (2006): 1562-1574.

22. H. Kvanme, Ø. Borgan, and I. Scheel, "Time-to-Event Prediction With Neural Networks and Cox Regression," *Journal of Machine Learning Research* 20 (2019): 1-30.

23. C. Lee, W. Zame, and J. Yoon, "DeepHit: A Deep Learning Approach to Survival Analysis With Competing Risks," 2018.

24. J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger, "DeepSurv: Personalized Treatment Recommender System Using a Cox Proportional Hazards Deep Neural Network," *BMC Medical Research Methodology* 18 (2018): 24.

25. M. F. Gensheimer and B. Narasimhan, "A Scalable Discrete-Time Survival Model for Neural Networks," *PeerJ* 7 (2019): e6257.

26. H. Kvanme and Ø. Borgan, "Continuous and Discrete-Time Survival Prediction With Neural Networks," *Lifetime Data Analysis* 27 (2021): 710-736.

27. O. Bouaziz and G. Nuel, "L0 Regularisation for the Estimation of Piecewise Constant Hazard Rates in Survival Analysis," (2016), https://arxiv.org/abs/1609.04595.

28. S. Pölsterl, N. Navab, and A. Katouzian, "Fast Training of Support Vector Machines for Survival Analysis," in *Machine Learning and Knowledge Discovery in Databases*, vol. 9285, ed. A. Appice (Cham, Switzerland: Springer International Publishing, 2015), 243-259.

29. S. Polsterl, N. Navab, and A. Katouzian, "An Efficient Training Algorithm for Kernel Survival Support Vector Machines," 2016.
## Comparison of different predictor models

The best model selected from the benchmarking step was used to assess the predictive quality of different set of predictors: baseline (BSL) features, BSL + BOR, BSL + TK and TK. These models were assessed on the training set using 10 repeats of 5-fold cross validations. The generalization quality of the models was assessed on the test set using 100 bootstrap sets. Furthermore, a Kruskal-Wallis test was used to compare the c-index of the models.

## Different Unbiased Two-step modeling approaches: TK4-OS4 and TKL1-PPS

To address the time-dependent covariate biases, one cannot rely on the TKL1 parameters to predict OS. Therefore, different unbiased two-step modeling approaches were investigated. On one hand, TK4-OS4, uses 4 cycles SLD data to estimate TK parameters (TK4) and predicts post cycle 4 OS. On the other hand, TKL1-PPS, uses TK parameters estimated from longitudinal SLD measurements restricted to first line + maintenance and predicts PPS.
30. T. Chen and C. Guestrin, "XGBoost: A Scalable Tree Boosting System," in *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining* (San Francisco, California, USA: Association for Computing Machinery, 2016), 785-794, https://doi.org/10.1145/2939672.2939785.

31. H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer, "Random Survival Forests," *Annals of Applied Statistics* 2 (2008): 841-860.

32. L. L. George, "The Statistical Analysis of Failure Time Data," *Technometrics* 45 (2003): 265-266.

33. R. Tibshirani, "The Lasso Method for Variable Selection in the Cox Model," *Statistics in Medicine* 16 (1997): 385-395.

34. G. Tutz and H. Binder, "Boosting Ridge Regression," *Computational Statistics & Data Analysis* 51 (2007): 6044-6059.

35. L. J. Wei, "The Accelerated Failure Time Model: A Useful Alternative to the Cox Regression Model in Survival Analysis," *Statistics in Medicine* 11 (1992): 1871-1879.

36. F. E. Harrell and R. M. Califf, "Evaluating the Yield of Medical Tests," *JAMA* 247 (1982): 2543-2546.

37. R. Bruno, D. Bottino, D. P. de Alwis, et al., "Progress and Opportunities to Advance Clinical Cancer Therapeutics Using Tumor Dynamic Models," *Clinical Cancer Research* 26 (2020): 1787-1795.

38. W. D. Stein, W. D. Figg, W. Dahut, et al., "Tumor Growth Rates Derived From Data for Patients in a Clinical Trial Correlate Strongly With Patient Survival: A Novel Strategy for Evaluation of Clinical Trial Data," *Oncologist* 13 (2008): 1046-1054.

39. L. Claret, P. Girard, J. O'Shaughnessy, et al., "Model-Based Predictions of Expected Anti-Tumor Response and Survival in Phase III Studies Based on Phase II Data of an Investigational Agent," *Journal of Clinical Oncology* 24 (2006): 6025.

40. Y. Wang, C. Sung, C. Dartois, et al., "Elucidation of Relationship Between Tumor Size and Survival in Non-Small-Cell Lung Cancer Patients Can Aid Early Decision Making in Clinical Drug Development," *Clinical Pharmacology and Therapeutics* 86 (2009): 167-174.

41. L. Claret, F. Mercier, B. E. Houk, P. A. Milligan, and R. Bruno, "Modeling and Simulations Relating Overall Survival to Tumor Growth Inhibition in Renal Cell Carcinoma Patients," *Cancer Chemotherapy and Pharmacology* 76 (2015): 567-573.

42. P. Chan, X. Zhou, N. Wang, Q. Liu, R. Bruno, and J. Y. Jin, "Application of Machine Learning for Tumor Growth Inhibition - Overall Survival Modeling Platform," *CPT: Pharmacometrics & Systems Pharmacology* 10 (2021): 59-66.

43. R. Bruno, M. Marchand, K. Yoshida, et al., "Tumor Dynamic Model-Based Decision Support for Phase Ib/II Combination Studies: A Retrospective Assessment Based on Resampling of the Phase III Study IMpower150," *Clinical Cancer Research* OF1-OF9 (2023): 1047-1055, https://doi.org/10.1158/1078-0432.CCR-22-2323.

44. J. G. Ibrahim, H. Chu, and L. M. Chen, "Basic Concepts and Methods for Joint Models of Longitudinal and Survival Data," *JCO* 28 (2010): 2796-2801.

45. A. Król, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels, and V. Rondeau, "Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event," *Journal of Statistical Software* 81 (2017): 1-52.

46. C. Tardivon, S. Desmée, M. Kerioui, et al., "Association Between Tumor Size Kinetics and Survival in Patients With Urothelial Carcinoma Treated With Atezolizumab: Implication for Patient Follow-Up," *Clinical Pharmacology and Therapeutics* 106 (2019): 810-820.

47. E. Christodoulou, J. Ma, G. S. Collins, E. W. Steyerberg, J. Y. Verbelak, and B. van Calster, "A Systematic Review Shows No Performance

Benefit of Machine Learning Over Logistic Regression for Clinical Prediction Models," *Journal of Clinical Epidemiology* 110 (2019): 12-22.

48. R. Knauer, M. Grimm, and E. Rodner, "PMLBmini: A Tabular Classification Benchmark Suite for Data-Scarce Applications," *arXiv* (2024), https://doi.org/10.48550/arXiv.2409.01635.

49. L. Grinsztajn, E. Oyallon, and G. Varoquaux, "Why Do Tree-Based Models Still Outperform Deep Learning on Typical Tabular Data?," in *Proceedings of the 36th International Conference on Neural Information Processing Systems* (Red Hook, NY, USA: Curran Associates Inc., 2024), 507-520.

50. E. Schindler, M. A. Amantea, M. O. Karlsson, and L. E. Friberg, "PK-PD Modeling of Individual Lesion FDG-PET Response to Predict Overall Survival in Patients With Sunitinib-Treated Gastrointestinal Stromal Tumor," *CPT: Pharmacometrics & Systems Pharmacology* 5 (2016): 173-181.

51. F. Mercier, M. Kerioui, S. Desmée, J. Guedj, O. Krieter, and R. Bruno, "Longitudinal Analysis of Organ-Specific Tumor Lesion Sizes in Metastatic Colorectal Cancer Patients Receiving First Line Standard Chemotherapy in Combination With Anti-Angiogenic Treatment," *Journal of Pharmacokinetics and Pharmacodynamics* 47 (2020): 613-625.

52. M. Kerioui, S. Desmée, F. Mercier, et al., "Assessing the Impact of Organ-Specific Lesion Dynamics on Survival in Patients With Recurrent Urothelial Carcinoma Treated With Atezolizumab or Chemotherapy," *ESMO Open* 7 (2022): 100346.

53. E. Schindler, S. M. Krishnan, R. Mathijssen, A. Ruggiero, G. Schia von, and L. E. Friberg, "Pharmacometric Modeling of Liver Metastases' Diameter, Volume, and Density and Their Relation to Clinical Outcome in Imatinib-Treated Patients With Gastrointestinal Stromal Tumors," *CPT: Pharmacometrics & Systems Pharmacology* 6 (2017): 449-457.

54. M. L. Maitland, J. Wilkerson, S. Karovic, et al., "Enhanced Detection of Treatment Effects on Metastatic Colorectal Cancer With Volumetric CT Measurements for Tumor Burden Growth Rate Evaluation," *Clinical Cancer Research* 26 (2020): 6464-6474.

55. J. Zhou, Q. Li, and Y. Cao, "Spatiotemporal Heterogeneity Across Metastases and Organ-Specific Response Informs Drug Efficacy and Patient Survival in Colorectal Cancer," *Cancer Research* 81 (2021): 2522-2533.

56. R. Almufti, M. Wilbaux, A. Oza, et al., "A Critical Review of the Analytical Approaches for Circulating Tumor Biomarker Kinetics During Treatment," *Annals of Oncology* 25 (2014): 41-56.

57. D. M. Kurtz, M. S. Esfahani, F. Scherer, et al., "Dynamic Risk Profiling Using Serial Tumor Biomarkers for Personalized Outcome Prediction," *Cell* 178 (2019): 699-713.

58. I. Netterberg, C. C. Li, L. Molinero, et al., "A PK/PDAnalysis of Circulating Biomarkers and Their Relationship to Tumor Response in Atezolizumab-Treated Non-Small Cell Lung Cancer Patients," *Clinical Pharmacology and Therapeutics* 105 (2018): 486-495.

59. I. Irurzun-Arana, E. Asín-Prieto, S. Martín-Algarra, and I. F. Trocóniz, "Predicting Circulating Biomarker Response and Its Impact on the Survival of Advanced Melanoma Patients Treated With Adjuvant Therapy," *Scientific Reports* 10 (2020): 7478.

60. K. H. Khan, D. Cunningham, B. Werner, et al., "Longitudinal Liquid Biopsy and Mathematical Modeling of Clonal Evolution Forecast Time to Treatment Failure in the PROSPECT-C Phase II Colorectal Cancer Clinical Trial," *Cancer Discovery* 8 (2018): 1270-1285.

61. L. N. Phuong, S. Salas, and S. Benzekry, "Computational Modeling Approaches for Circulating Cell-Free DNA in Oncology," 2024.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section.
# RESULTS

## Patients and disease characteristics

Overall, 526 patients treated with either the EXTREME or TPEx regimen were retrieved from the GORTEC 2014-01 TPExtreme clinical trial⁵. Patients and disease characteristics are summarized in Table S3. 60.1% of the patients received a second line treatment while 33.2% got a third line treatment and finally 13.3% received a fourth line treatment. 44.6% of the primary tumors were in the oropharynx and 64% of the patients had metastasis at inclusion (Table S3). The typical treatment history of a patient is depicted in **Figure 1.A**. It highlights the inter-patient heterogeneity in terms of both efficacy and toxicity, thus emphasizing the complexity of clinical management of such patients at bedside.

## Tumor kinetics modeling

The SLD measurements on 368 (70%) patients – restricted to 4 cycles (total number of time points = 672 and median number of time points per patient = 2) – were used to calibrate the different TK models. All tested models showed a good agreement with the data. However, according to the corrected Bayesian Information Criterion (BICc), model 1 was found to be the best to describe the data in the different cases (TK4 and TKL1) (Table S1 for TK4). All the parameters were identifiable with a maximum RSE of 8.73% (Table S2).

Diagnostic plots did not reveal any misspecification of the model at the individual level (Figure 2. A-B). Recalibration of the model using the treatment arm as a covariate to the kinetic parameters (KS, KG) revealed a discernible effect of the TPEx treatment on the tumor growth rate (KG) in comparison to its effect on the tumor shrinkage rate (KS). A Wald test demonstrated a statistically significant association between the treatment arm and KG, ($p = 6.35 \times 10^{-3}$). No association between the treatment arm and KS was found. Additionally, simulations conducted using the estimated population parameters suggested a slightly longer mean time to progression in the TPEx group compared to the EXTREME group (Figure 2. C). However, adding any covariate (including treatment arm) to the model resulted in a decreased goodness-of-fit and triggered identifiability issues. Therefore, a model without covariates was kept for the predictive analysis, on the entire patient cohort (i.e., both treatment arms).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Diagnostic plots and treatment arms simulation. (A) Representative individual fits. The purple line represents the model prediction and the blue dots the observed SLDs. The red, green and dark-blue dotted vertical lines represent the end of the first, second line and third lines, respectively. (B) Observed SLDs vs predicted SLDs (C) Simulations of tumor kinetics of each treatment arm using the typical population parameters (fixed-effects). SLD: sum of largest diameters</figcaption></figure>

ML for survival analysis

ML models benchmarking: unbiased predictions using TK4–OS4 et TK1L–PPS.

To avoid time-dependent covariate bias, we considered either predicting post-cycle 4 survival (OS4) from data restricted to the first 4 treatment cycles (TK4–OS4) or predict post-progression survival (PPS) from TK parameters derived from the data during the first line of treatment until progression (TKL1–PPS). On the training set, for both TK4–OS4 – and TKL1–PPS, the random survival forest (RSF) model did not significantly outperformed Cox PH or AFT (median c-indices = 0.65, 0.635 and 0.63 for TK4–OS4, **Figure 3** and 0.601, 0.598 and 0.595 for TKL1–PPS, Figure S1). However, an examination of the results on the test set for TK4–OS4 revealed that RSF exhibited the best generalization properties in comparison to Cox PH and AFT (**Figure S2**). Overall, TKL1-PPS (best training c-index = 0.6) exhibited worse performances than TK4–OS4 (best training c-index = 0.65, Figure 3B). Moreover, deep
learning architectures, specifically DeepHit and DeepSurv, exhibited suboptimal performances for both TK4-OS4 and TKL1-PPS. Upon examination of the test set, for TK4-OS4, the Cox PH model achieved a mean c-index of 0.54 (95% CI: 0.52 - 0.56) **(Figure S2)**.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Benchmark of ML methods, TK4 – OS4. (A) Graphical representation of the different ML and statistical methods. The green boxes are the tested models. (B) Boxplots of the results of the benchmarking process on OS4 prediction. The red Diamond is the median c-index and the middle line the mean c-index.</figcaption></figure>

## Feature importance.

Noteworthy, the most important feature for OS4 prediction was first the patient performance status followed by the TK parameters (*logKG*, *logKS* and *TTG*) and the baseline tumor size *TS*₀ (**Figure 4. A**). In addition, utilizing the maximally selected rank statistic, *TS*₀, *logKG* and *logKS* provided a significant separation of the survival curves (log-rank *p* < 0.01 and *p* < 0.05 respectively) (**Figure 4. B-D**).

However, for TKL1-PPS, the top-3 parameters were loco-regional relapse, followed by *TS*₀, the performance status (**Figure S3.A**). Surprisingly, Unlike TK4 – OS4, *logKS*, triggered crossed survival curves (**Figure S3.D**) potentially indicating a time-dependent nature of the growth rate parameter.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Features Importance and stratified survival curves for TK4-OS4. (A) Feature importance based on the random survival forest model and a permutation importance algorithm. (B) Survival curves with optimal cut-off on *TS*₀ (SLD at diagnosis) using the maximally selected rank statistics from the 'maxstat' R package. (C) Survival curves with optimal cut-off on *logKG* (tumor growth rate). (D) Survival curves with optimal cut-off on *logKS* (tumor shrinkage rate).</figcaption></figure>

Comparison of different predictor models: TK4 parameters are better predictors of OS4 than BOR.

To compare the relative predictive performances of BOR and TK4 as predictors of OS4, four sets of predictors were considered in combination with RSF: 12 baseline features (BSL, see **Table S1**) alone, BSL + BOR, BSL + TK4 and TK4 (Figure 5). In the training set, the BSL + TK4 outperformed BOR + BSL (0.63 (95% CI: 0.61 – 0.65) versus 0.612
(95% CI: 0.58 - 0.64), respectively, *p* < 0.001, Kruskal-Wallis, **Figure 5.A**). This trend persisted in the test dataset, with respective c-index of 0.63 (95% CI: 0.61 - 0.65) and 0.52 (95% CI: 0.50 - 0.54) (*p* < 0.0001, **Figure 5.B**). Intriguingly, the integration of BSL with BOR adversely impacted the model accuracy on test set as compared to BSL alone (**Figure 5.B**). However, this was not observed in TKL1-PPS, such behaviour was absent (**Figure S4.A-B**).

<figure><img src="image_7.png" /><figcaption>Figure 5. Comparison of different predictive models for TK4-OS4.</figcaption></figure>

The red diamond and middle line are the median and mean c-index, respectively. (A) Predictive models on the training set. (B) Predictive models on the test set. *p* refers to the p-value of the Kruskal-Wallis test.
# DISCUSSION

Tumor kinetics modeling to predict overall survival (TK–OS) has been extensively developed in the last 15 years³². Multiple cancer types, often in the advanced or metastatic stage have been studied, including prostate³³, breast³⁴, colorectal⁹, non-small cell lung cancer³⁵ or renal cell carcinoma³⁶. Our study is the first to apply such modeling to HNSCC, using data from the TPExtreme trial⁵. In addition, it is one of the few to consider and benchmark machine learning survival models versus classical (semi-)parametric survival analysis methods³⁷.

As often observed for the other cancer types, on-treatment TK could be accurately described by a simple pattern mathematically: the sum of two exponentials. A decreasing phase governed by a shrinkage parameter $KS_{g}$ and a regrowth phase led by a parameter $KG$. Using the treatment arm as a covariate showed a significantly slower regrowth (parameter $KG$) and deeper SLD decrease in the investigational (TPEx) arm versus the control (EXTREME). However, even though the decrease in SLD is deeper in the TPEx arm, this decrease is not necessarily faster in this arm compared to the EXTREME arm, at least within the first 3 months.

This result thus prompts discussion on the necessity of undergoing TPEx for patients who present with a high tumor burden or symptoms that require a rapid response⁵. Additionally, observations of simulated median TK suggested a benefit of TPEx versus EXTREME, observed at ~ 12 months, with a better safety profile⁵.

Unlike several previous TK modeling studies mostly aimed to assist drug development⁹,¹²,³⁴,³⁵,³⁸, here we turned our interest on the ability of TK metrics to predict earlier individual survival. In such context, the value of predictive modeling does not rely on the ability to forecast study-level survival curves but rather lies in the evaluation of discrimination metrics such as the c-index. It is also crucial to avoid time-dependent covariate bias³⁹,⁴⁰. One approach to do so is to rely on joint modeling⁴¹⁻⁴³. Here, motivated by clinically relevant scenarios in routine care, we rather focused on two operational survival prediction problems. The first was to evaluate the survival predictive value of on-treatment TK metrics during the four cycles corresponding to TPEx chemotherapy duration (TK4–OS4). The second leveraged data from the first line of treatment plus maintenance to predict post-progression survival (TK1L–PPS). These could be of value to guide clinical decisions regarding maintenance or second line therapy. The latter setting is of particular interest because, while clear guidelines are established from large phase 3 trials in the first line setting, the diversity of subsequent lines options is wider and often left to the appreciation of the medical oncologist. Disposing of an individualized PPS predictive tool could for example help decide whether to orientate a given patient to one therapeutic option or another and thus optimize the treatment sequence.
To maximize predictive performances and account for possible nonlinear interactions between
predictors of survival, we benchmarked the current state-of-the art survival machine learning
models in our setting. We therefore compared traditional (Cox proportional hazard) modeling
to tree ensemble models (e.g., gradient boosting or random survival forest) and neural
networks (CoxTime, DeepSurv, DeepHit). Our results indicate that the more complex (and less
interpretable) models only rarely and marginally surpassed the linear standard (Cox) model.
These findings align with other studies showing that complex ML might not overcome classical
statistical models for cases relative to predictive modeling in healthcare from structured data
⁴⁴. Rather, advanced neural network engineering performs better for computer vision or natural
language processing tasks.

Yet, the best model for TK4-OS4 was a random survival forest model. It demonstrated a
modest (but significant) predictive performance of individual survival with a c-index of 0.592.
When reducing the number of variables to only BSL + TK4, the model: 1) did not exhibit signs
of overfitting (same performances in the training and test sets) but most importantly 2)
indicated a substantially better prediction of OS4 than a model using BSL + BOR. This last
point indicates that the current standard RECIST guidelines for evaluation and monitoring of
response could be improved using dynamic TK-based metrics that account for the entire on-
treatment longitudinal data. Indeed, RECIST is only based on the mere static difference
between the current time point and the baseline or nadir.

Unfortunately, our results indicated a limited value of TKL1 model-based metrics for PPS
prediction. However, several refinements could be further investigated to improve. For
example, modeling individual lesion sizes rather than the SLD unraveled substantial inter-
lesion and inter-organ variability within a specific patient, and that this variability was linked to
response ⁴⁵⁻⁴⁷. In addition, considering the total volume rather one-dimensional measurements
can reveal inter-lesion variability that could not be observed with diameters ⁴⁸ and provide
better assessment of the treatment effect ⁴⁹. Spatiotemporal heterogeneity across metastases
has also been found relevant in a large study comprising more than 11,000 lesions in 2,802
colorectal cancer patients ⁵⁰. Further, modeling the kinetics of functional tumor markers could
add predictive value. For instance, Schindler et al. used the per-lesion maximal standard
uptake value (SUVmax) obtained from positron emission tomography and not only found
significant inter-lesion variability but also that this marker outperformed tumor size for
prediction of survival ⁴⁵.

Beyond invasive or irradiating tumor markers, circulating markers could also enhance
predictions ⁵¹⁻⁵⁴. In a study modeling 4 simple blood markers kinetics (BK) from routine lab
tests in addition to TK in a large dataset (862 patients in the training set for model development,
553 in the test set for model evaluation), and combining mixed-effects modeling with ML, we
recently found significant added OS predictive value of BK parameters compared with baseline
+ TK parameters <sup>16,18</sup>. One of the main advantages of such markers is that they are non-
invasive and can be sampled more frequently than tumor measurements. In addition, the
advent of liquid biopsies provides qualitative and quantitative data on circulating free or tumor
DNA that pave a promising avenue to develop more accurate predictive computational models
leveraging on-treatment data <sup>55,56</sup>.

In conclusion, our work suggests that early analysis of first-line tumor kinetics by ML modeling,
likely reflecting tumor biology, could eventually enable the individual assessment of post-
treatment survival. This could also serve as a tool to guide the selection of different therapeutic
options and thus aid in studying various treatment sequences.
## ACKNOWLEDGEMENTS

This work is part of the QUANTIC project funded by ITMO Cancer AVIESAN and the French Institute National du Cancer (grant #19CM148-00)
REFERENCES

1. Johnson, D. E. et al. Head and neck squamous cell carcinoma. *Nat Rev Dis Primers* **6**, 92 (2020).

2. Okano, S. et al. Induction chemotherapy in locally advanced squamous cell carcinoma of the head and neck. *Japanese Journal of Clinical Oncology* **51**, 173–179 (2021).

3. Vermorken, J. B. et al. Platinum-Based Chemotherapy plus Cetuximab in Head and Neck Cancer. *N Engl J Med* **359**, 1116–1127 (2008).

4. Burtness, B. et al. Pembrolizumab alone or with chemotherapy versus cetuximab with chemotherapy for recurrent or metastatic squamous cell carcinoma of the head and neck (KEYNOTE-048): a randomised, open-label, phase 3 study. *The Lancet* **394**, 1915–1928 (2019).

5. Guigay, J. et al. Cetuximab, docetaxel, and cisplatin versus platinum, fluorouracil, and cetuximab as first-line treatment in patients with recurrent or metastatic head and neck squamous-cell carcinoma (GORTEC 2014-01 TPExtreme): a multicentre, open-label, randomised, phase 2 trial. *The Lancet Oncology* **22**, 463–475 (2021).

6. Eisenhauer, E. A. et al. New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1). *European Journal of Cancer* **45**, 228–247 (2009).

7. Saad, E., Katz, A., Machado, K. & Buyse, M. Post-Progression Survival (PPS) and Overall Survival (OS) According to Treatment Type in Contemporary Phase III Trials in Advanced Breast Cancer (ABC). *Cancer Research* **69**, 5116–5116 (2009).

8. Imai, H. et al. Progression-free survival, post-progression survival, and tumor response as surrogate markers for overall survival in patients with extensive small cell lung cancer. *Annals of Thoracic Medicine* **10**, (2015).

9. Claret, L. et al. Model-Based Prediction of Phase III Overall Survival in Colorectal Cancer on the Basis of Phase II Tumor Dynamics. *Journal of Clinical Oncology* **27**, 4103–4108 (2009).
10. Claret, L. et al. Evaluation of Tumor-Size Response Metrics to Predict Overall Survival in Western and Chinese Patients With First-Line Metastatic Colorectal Cancer. *Journal of Clinical Oncology* **31**, 2110–2114 (2013).

11. Claret, L. et al. Evaluation of Tumor-Size Response Metrics to Predict Overall Survival in Western and Chinese Patients With First-Line Metastatic Colorectal Cancer. *Journal of Clinical Oncology* **31**, 2110–2114 (2013).

12. Claret, L. et al. A Model of Overall Survival Predicts Treatment Outcomes with Atezolizumab versus Chemotherapy in Non-Small Cell Lung Cancer Based on Early Tumor Kinetics. *Clinical Cancer Research* **24**, 3292–3298 (2018).

13. Beyersmann, J., Gastmeier, P., Wolkewitz, M. & Schumacher, M. An easy mathematical proof showed that time-dependent bias inevitably leads to biased effect estimation. *J Clin Epidemiol* **61**, 1216–1221 (2008).

14. Beyersmann, J., Wolkewitz, M. & Schumacher, M. The impact of time-dependent bias in proportional hazards modelling. *Stat Med* **27**, 6439–6454 (2008).

15. Mistry, H. On the relationship between tumour growth rate and survival in non-small cell lung cancer. (2017).

16. Benzekry, S. et al. Predicting survival and trial outcome in non-small cell lung cancer integrating tumor and blood markers kinetics with machine learning. 2023.09.26.23296135 Preprint at https://doi.org/10.1101/2023.09.26.23296135 (2024).

17. Benzekry, S. Artificial Intelligence and Mechanistic Modeling for Clinical Decision Making in Oncology. *Clin. Pharmacol. Ther.* **108**, 471–486 (2020).

18. Benzekry, S. et al. Supporting decision making and early prediction of survival for oncology drug development using a pharmacometrics-machine learning based model. in *PAGE* vol. 30 10276 (2022).

19. Stein, W. D., Yang, J., Bates, S. E. & Fojo, T. Bevacizumab Reduces the Growth Rate Constants of Renal Carcinomas: A Novel Algorithm Suggests Early Discontinuation of Bevacizumab Resulted in a Lack of Survival Advantage. *The Oncologist* **13**, 1055–1062 (2008).
20. Samson, A., Lavielle, M. & Mentré, F. Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model. *Computational Statistics & Data Analysis* **51**, 1562–1574 (2006).

21. Kvamme, H., Borgan, Ø. & Scheel, I. Time-to-Event Prediction with Neural Networks and Cox Regression. Preprint at (2019).

22. Lee, C., Zame, W. & Yoon, J. DeepHit: A Deep Learning Approach to Survival Analysis with Competing Risks. **8**.

23. Katzman, J. L. et al. DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network. *BMC Med Res Methodol* **18**, 24 (2018).

24. Gensheimer, M. F. & Narasimhan, B. A scalable discrete-time survival model for neural networks. *PeerJ* **7**, e6257 (2019).

25. Kvamme, H. & Borgan, Ø. Continuous and discrete-time survival prediction with neural networks. *Lifetime Data Analysis* **27**, 710–736 (2021).

26. Pölsterl, S., Navab, N. & Katouzian, A. Fast Training of Support Vector Machines for Survival Analysis. in *Machine Learning and Knowledge Discovery in Databases* (eds. Appice, A. et al.) vol. 9285 243–259 (Springer International Publishing, Cham, 2015).

27. Polsterl, S., Navab, N. & Katouzian, A. An Efficient Training Algorithm for Kernel Survival Support Vector Machines. **13**.

28. Ishwaran, H., Kogalur, U. B., Blackstone, E. H. & Lauer, M. S. Random survival forests. *The Annals of Applied Statistics* **2**, (2008).

29. Tutz, G. & Binder, H. Boosting ridge regression. *Computational Statistics & Data Analysis* **51**, 6044–6059 (2007).

30. Wei, L. J. The accelerated failure time model: A useful alternative to the cox regression model in survival analysis. *Statistics in Medicine* **11**, 1871–1879 (1992).

31. Harrell, F. E. & Califf, R. M. Evaluating the Yield of Medical Tests. **4**.

32. Bruno, R. et al. Progress and Opportunities to Advance Clinical Cancer Therapeutics Using Tumor Dynamic Models. *Clinical Cancer Research* **26**, 1787–1795 (2020).
33. Stein, W. D. et al. Tumor growth rates derived from data for patients in a clinical trial correlate strongly with patient survival: a novel strategy for evaluation of clinical trial data. *The Oncologist* **13**, 1046–1054 (2008).

34. Claret, L. et al. Model-based predictions of expected anti-tumor response and survival in phase III studies based on phase II data of an investigational agent. *J Clin Oncol* (2006) doi:10.1200/jco.2006.24.18_suppl.6025;issue:issue:10.1200/jco.2006.24 issue-18_suppl.

35. Wang, Y. et al. Elucidation of relationship between tumor size and survival in non-small-cell lung cancer patients can aid early decision making in clinical drug development. *Clin Pharmacol Ther* **86**, 167–174 (2009).

36. Claret, L., Mercier, F., Houk, B. E., Milligan, P. A. & Bruno, R. Modeling and simulations relating overall survival to tumor growth inhibition in renal cell carcinoma patients. *Cancer Chemother Pharmacol* **76**, 567–573 (2015).

37. Chan, P. et al. Application of Machine Learning for Tumor Growth Inhibition - Overall Survival Modeling Platform. *CPT Pharmacometrics Syst Pharmacol* **10**, 59–66 (2021).

38. Bruno, R. et al. Tumor Dynamic Model-Based Decision Support for Phase Ib/II Combination Studies: A Retrospective Assessment Based on Resampling of the Phase III Study IMpower150. *Clinical Cancer Research* OF1–OF9 (2023) doi:10.1158/1078-0432.CCR-22-2323.

39. Mistry, H. Time-Dependent bias of tumor growth rate and time to tumor regrowth. *CPT: Pharmacometrics & Systems Pharmacology* **5**, 587–587 (2016).

40. van Walraven, C., Davis, D., Forster, A. J. & Wells, G. A. Time-dependent bias was common in survival analyses published in leading clinical journals. *J Clin Epidemiol* **57**, 672–682 (2004).

41. Ibrahim, J. G., Chu, H. & Chen, L. M. Basic Concepts and Methods for Joint Models of Longitudinal and Survival Data. *JCO* **28**, 2796–2801 (2010).

42. Król, A. et al. Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event. *Journal of Statistical Software* **81**, 1–52 (2017).
43. Tardivon, C. et al. Association Between Tumor Size Kinetics and Survival in Patients With Urothelial Carcinoma Treated With Atezolizumab: Implication for Patient Follow-Up. *Clin Pharmacol Ther* **106**, 810–820 (2019).

44. Christodoulou, E. et al. A systematic review shows no performance benefit of machine learning over logistic regression for clinical prediction models. *Journal of Clinical Epidemiology* **110**, 12–22 (2019).

45. Schindler, E., Amantea, M. A., Karlsson, M. O. & Friberg, L. E. PK-PD modeling of individual lesion FDG-PET response to predict overall survival in patients with sunitinib-treated gastrointestinal stromal tumor. *CPT: Pharmacometrics & Systems Pharmacology* **5**, 173–181 (2016).

46. Mercier, F. et al. Longitudinal analysis of organ-specific tumor lesion sizes in metastatic colorectal cancer patients receiving first line standard chemotherapy in combination with anti-angiogenic treatment. *J Pharmacokinet Pharmacodyn* **47**, 613–625 (2020).

47. Kerioui, M. et al. Assessing the impact of organ-specific lesion dynamics on survival in patients with recurrent urothelial carcinoma treated with atezolizumab or chemotherapy. *ESMO Open* **7**, 100346 (2022).

48. Schindler, E. et al. Pharmacometric Modeling of Liver Metastases’ Diameter, Volume, and Density and Their Relation to Clinical Outcome in Imatinib-Treated Patients With Gastrointestinal Stromal Tumors. *CPT: Pharmacometrics & Systems Pharmacology* **6**, 449–457 (2017).

49. Maitland, M. L. et al. Enhanced Detection of Treatment Effects on Metastatic Colorectal Cancer with Volumetric CT Measurements for Tumor Burden Growth Rate Evaluation. *Clinical Cancer Research* **26**, 6464–6474 (2020).

50. Zhou, J., Li, Q. & Cao, Y. Spatiotemporal Heterogeneity across Metastases and Organ-Specific Response Informs Drug Efficacy and Patient Survival in Colorectal Cancer. *Cancer Research* **81**, 2522–2533 (2021).

51. Almufti, R. et al. A critical review of the analytical approaches for circulating tumor biomarker kinetics during treatment. *Ann Oncol* **25**, 41–56 (2014).
52. Kurtz, D. M. et al. Dynamic Risk Profiling Using Serial Tumor Biomarkers for Personalized Outcome Prediction. *Cell* **178**, 699-713.e19 (2019).

53. Netterberg, I. et al. A PK/PDAnalysis of Circulating Biomarkers and Their Relationship to Tumor Response in Atezolizumab-Treated non-small Cell Lung Cancer Patients. *Clin Pharmacol Ther* **105**, 486-495 (2018).

54. Irurzun-Arana, I., Asín-Prieto, E., Martín-Algarra, S. & Trocóniz, I. F. Predicting circulating biomarker response and its impact on the survival of advanced melanoma patients treated with adjuvant therapy. *Sci Rep* **10**, 7478 (2020).

55. Khan, K. H. et al. Longitudinal Liquid Biopsy and Mathematical Modeling of Clonal Evolution Forecast Time to Treatment Failure in the PROSPECT-C Phase II Colorectal Cancer Clinical Trial. *Cancer Discovery* **8**, 1270-1285 (2018).

56. Phuong, L. N., Salas, S. & Benzekry, S. Computational modeling approaches for circulating cell-free DNA in oncology. (2024).
# SUPPLEMENTARY

<table><caption>Table S1 Goodness-of-fit of the different TGI models on TK4</caption><thead><tr><th>Models</th><th>AIC</th><th>BIC</th><th>BICc</th><th>OFV</th></tr></thead><tbody><tr><td>Model 1</td><td><strong>11896.57</strong></td><td><strong>11916.11</strong></td><td><strong>11920.61</strong></td><td><strong>11886.57</strong></td></tr><tr><td>Model 2</td><td>12337.5</td><td>12364.86</td><td>12370.86</td><td>12323.5</td></tr><tr><td>Model 3</td><td>11919.28</td><td>11946.64</td><td>11952.64</td><td>11905.28</td></tr></tbody></table>

<table><caption>Table S2 Population parameter estimates</caption><thead><tr><th>Parameter</th><th>Value</th><th>S.E.</th><th>R.S.E. (%)</th></tr></thead><tbody><tr><td>KS<sub>pop</sub></td><td>0.27</td><td>0.0163</td><td>6.03</td></tr><tr><td>KG<sub>pop</sub></td><td>0.0309</td><td>0.0027</td><td>8.73</td></tr><tr><td>ω<sub>KS</sub></td><td>0.821</td><td>0.0544</td><td>6.63</td></tr><tr><td>ω<sub>KG</sub></td><td>1</td><td>0.0629</td><td>6.26</td></tr><tr><td>a</td><td>7.43</td><td>0.303</td><td>4.08</td></tr></tbody></table>

<table><caption>Table S3 Patients and disease characteristics</caption><thead><tr><th></th><th>n (%) or mean ± SD (min-max)</th></tr></thead><tbody><tr><td><strong>Sex</strong></td><td></td></tr><tr><td>Female</td><td>66 (12.54%)</td></tr><tr><td>Male</td><td>460 (87.46%)</td></tr><tr><td><strong>Age (years)</strong></td><td>59.6 ± 6.5 (23.3 – 71.4)</td></tr><tr><td><strong>Alcohol consumption</strong></td><td></td></tr><tr><td>No</td><td>133 (25.29%)</td></tr><tr><td>Yes (past)</td><td>206 (39.16%)</td></tr><tr><td>Yes (currently)</td><td>182 (34.6%)</td></tr><tr><td><strong>Tobacco status</strong></td><td></td></tr><tr><td>No</td><td>41 (7.8%)</td></tr><tr><td>Yes (past)</td><td>327 (62.17%)</td></tr><tr><td>Yes (currently)</td><td>158 (30.04%)</td></tr><tr><td><strong>Primary tumor location</strong></td><td></td></tr><tr><td>Oropharynx</td><td>235 (44.68%)</td></tr><tr><td>Others</td><td>291 (55.32%)</td></tr><tr><td><strong>Disease evolution at inclusion</strong></td><td></td></tr><tr><td>Metastasis</td><td>337 (64.07%)</td></tr><tr><td>Locoregional relapse</td><td>305 (57.98%)</td></tr><tr><td>Metastasis and locoregional relapse</td><td>116 (22.05%)</td></tr></tbody></table>
**Previous cancer treatment**

<table><tr><td>No</td><td>90 (17.12%)</td></tr><tr><td>Yes</td><td>435 (82.7%)</td></tr></table>

**Further treatment lines**

<table><tr><td>Line 2</td><td>316 (60.07 %)</td></tr><tr><td>Line 3</td><td>175 (33.27 %)</td></tr><tr><td>Line 4</td><td>70 (13.31 %)</td></tr></table>

**ECOG Performance status**

<table><tr><td>1</td><td>351 (66.7 %)</td></tr><tr><td>0</td><td>175 (33.3 %)</td></tr></table>

**lymphadenectomy for residual lymph node**

<table><tr><td>yes</td><td>378 (71.8 %)</td></tr><tr><td>no</td><td>146 (27.8 %)</td></tr></table>

**TNM classification at inclusion**

T, size of the primary tumor

<table><tr><td>T0</td><td>2 (0.380 %)</td></tr><tr><td>T1</td><td>56 (10.6 %)</td></tr><tr><td>T2</td><td>154 (29.3 %)</td></tr><tr><td>T3</td><td>143 (27.2 %)</td></tr><tr><td>T4</td><td>151 (28.7 %)</td></tr><tr><td>Tx</td><td>14 (2.66 %)</td></tr></table>

N, extent to regional lymph nodes

<table><tr><td>N0</td><td>129 (24.5 %)</td></tr><tr><td>N1</td><td>69 (13.1 %)</td></tr><tr><td>N2</td><td>260 (49.4 %)</td></tr><tr><td>N3</td><td>48 (9.13 %)</td></tr><tr><td>Nx</td><td>14 (2.66 %)</td></tr></table>

M, presence of distant metastasis

<table><tr><td>M0</td><td>395 (75.1 %)</td></tr><tr><td>M1</td><td>96 (18.3 %)</td></tr><tr><td>Mx</td><td>28 (5.32 %)</td></tr></table>

**Squamous cell carcinoma (SCC) histology**

<table><tr><td>moderately differentiated</td><td>211 (40.1%)</td></tr><tr><td>poorly or undifferentiated</td><td>106 (20.2 %)</td></tr><tr><td>Well-differentiated</td><td>154 (29.3 %)</td></tr></table>
<table><caption>Table S4 BSL and TK features</caption><tr><td>Baseline features (BSL)</td><td><ul><li>TS<sub>0</sub></li><li>Age</li><li>Sex</li><li>Smoking status</li><li>Alcohol consumption</li><li>Primary tumor location</li><li>Metastases relapse</li><li>TNM classification at inclusion</li><li>SCC histology</li><li>Loco-regional relapse</li><li>Lymphadenectomy of the residual lymph node</li><li>Performance status</li></ul></td></tr><tr><td>Tumor Kinetic parameters (TK)</td><td><ul><li>logKG</li><li>logKS</li><li>TTG</li></ul></td></tr></table>

<figure><img src="image_8.png" /><figcaption>Figure S1 Benchmark of ML methods on the training set, TKL1-PPS model. Boxplots of the results of the benchmarking process on PPS prediction. The red Diamond is the median c-index and the middle line the mean c-index.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure S2 Best ML models on the test set, TK4 - OS4</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure S3 Parameters Importance and survival curves, TKL1-PPS model. (A) Feature importance based on the random survival forest model and a permutation importance algorithm. (B) Survival curves with optimal cut-off on TS₀ (SLD at diagnosis) using the maximally selected rank statistics from the 'maxstat' R package. (C) Survival curves with optimal cut-off on KG (tumor growth rate). (D) Survival curves with optimal cut-off on KS (tumor shrinkage rate)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure S4 Comparison of different predictors models, TKL1-PPS model.** The red Diamond is the median c-index and the middle line the mean c-index. *p* refers to the p-value of the Kruskal-Wallis test (A) Boxplots of the results of the comparison of the predictor models on the train set. (B) Boxplots of the results of the comparison of the predictor models on the test set.</figcaption></figure>