Predicting survival and trial outcome
in non-small cell lung cancer
integrating tumor and blood
markers kinetics with machine
learning

Sébastien Benzekry<sup>1</sup> <font color="red">1</font>✉, Mélanie Karlsen<sup>1</sup>, Célestin Bigarré<sup>1</sup>, Abdessamad El Kaoutari<sup>1</sup>, Bruno Gomes<sup>2</sup>, Martin Stern<sup>3</sup>, Ales Neubert<sup>4</sup>, Rene Bruno<sup>5</sup>, François Mercier<sup>6</sup>, Suresh Vatakuti<sup>7</sup>, Peter Curle<sup>8</sup>, Candice Jamois<sup>9</sup>

¹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; ²Pharma Research and Early Development, Early Development Oncology, Roche Innovation Center Basel, Switzerland; ³Pharma Research and Early Development, Early Development Oncology, Roche Innovation Center Zurich, Switzerland; ⁴Pharma Research and Early Development, Data & Analytics, Roche Innovation Center Basel, Switzerland; ⁵Modeling and Simulation, Clinical Pharmacology, Genentech Research and Early Development, Marseille France; ⁶Modeling and Simulation, Clinical Pharmacology, Genentech Research and Early Development, Roche Innovation Center Basel; ⁷Pharma Research and Early Development, Predictive Modeling and Data Analytics, Roche Innovation Center Basel, Switzerland; ⁸Invigate, Basel, Switzerland; ⁹Pharma Research and Early Development, Translational PKPD and Clinical Pharmacology, Roche Innovation Center Basel, Switzerland

Abstract

Existing survival prediction models rely only on baseline or tumor kinetics data and lack machine
learning integration. We introduce a novel kinetics-machine learning (kML) model that integrates
baseline markers, tumor kinetics and four on-treatment simple blood markers (albumin, CRP,
lactate dehydrogenase and neutrophils). Developed for immune-checkpoint inhibition (ICI) in
non-small cell lung cancer on three phase 2 trials (533 patients), kML was validated on the two
arms of a phase 3 trial (ICI and chemotherapy, 377 and 354 patients). It outperformed the current
state-of-the-art for individual predictions with a test set c-index of 0.790, a 12-months survival
accuracy of 78.7% and a hazard ratio of 25.2 (95% CI: 10.4 – 61.3, p < 0.0001) to identify long-term
survivors. Critically, kML predicted the success of the phase 3 trial using only 25 weeks of
on-study data (predicted HR = 0.814 (0.64 – 0.994) versus final study HR = 0.778 (0.65 – 0.931)). Our
model constitutes a valuable approach to support personalized medicine and drug development.

☑ For correspondence:
sebastien.benzekry@inria.fr

**Present address:** COMPO
team, Pharmacy faculty
27 Bd Jean Moulin
13385 Marseille, FRANCE

**Data availability:** Qualified researchers may request access to individual patient level data through the clinical study data request platform (https://vivli.org/). Further details on Roche's criteria for eligible studies are available here (https://vivli.org/members/ourmembers/). For further details on Roche's Global Policy on the Sharing of Clinical Information and how to request access to related clinical study documents, see here https://www.roche.com/innovation/process/clinical-trials/data-sharing/.

**Funding:** This work was sponsored by the Roche Pharma Research and Early Development (pRED) One-D Modeling and Simulation Digital Initiative. It also benefited from funding from ITMO Cancer AVIESAN and French Institut National du Cancer (grant #19CM148-00)

**Competing interests:** The authors declare the existence of a financial competing interest
## Introduction

Lung cancer is the leading cause of cancer death worldwide¹, with non-small cell lung cancer (NSCLC) being the most prevalent type, representing 80%–85% of case². Immune-checkpoint inhibitors (ICI) (e.g., atezolizumab (ATZ)) have led to significant improvements in survival rates for patients with advanced cancers such as NSCLC³,⁴. However, there is still a large variability in clinical response and progression eventually occurs in a majority of patients⁵. Additionally, drug development in immuno-oncology is highly challenging, with a 95% attrition rate⁶. Current approaches for go/no-go decisions are based on interim endpoints (e.g., progression-free survival, overall response rate) that have often been found to be poor predictors of the primary endpoint of most clinical trials in oncology, overall survival (OS)⁷. This calls for better surrogate markers at interim analyses. Altogether, there is a need for better and validated predictive models of OS for both personalized health care (individual predictions) and drug development (trial predictions).

Currently, PDL1 expression is the only routine biomarker used for NSCLC patients⁵,⁸ despite being controversial¹,¹⁰. Tumor mutational burden⁸,¹¹,¹² and transcriptomic data⁵,¹³,¹⁴ have also been investigated but did not reach clinical practice. Here we posit that such static and single marker approach is intrinsically limited and that substantial additional predictive performances could be gained by: 1) using multi-modal integrative analyses relying on a combination of markers and machine learning algorithms⁵,¹²,¹⁴,¹⁵ and 2) including dynamic markers obtained from early on-treatment data¹⁵,¹⁶. The nonlinear mixed-effects (NLME) modeling approach is well suited for the latter¹⁷, and tumor kinetics (TK) model-based metrics have been shown to carry significant predictive value for OS in oncology, including ATZ monotherapy in advanced NSCLC¹⁸⁻²⁰. The first main novelty of the current study is to establish the predictive value of model-based parameters of simple blood markers kinetics (BK), in addition to TK.

The second main novelty is to apply machine learning (ML) algorithms, increasingly used in biology and medicine²¹ but only rarely for TK-OS modeling²², instead of classical survival models. Extensions of classical ML models to survival data have been proposed (e.g., random survival forests²³), but their actual superiority over standard approaches remains controversial²⁴. In addition, most ML studies to date are underpowered due to low sample sizes in both training and test sets.

Here, we coupled the strengths of NLME modeling with ML to derive a predictive model of OS from baseline and on-treatment data, called kinetics-machine learning (kML, Figure 1A). We leveraged large training and test datasets to achieve robust results (Figure 1B). Subsequently, we tested the operational predictive capabilities of kML in two relevant scenarios: 1) individual prediction of OS and 2) prediction of the outcome of a phase 3 trial from early on-study data.

## Methods

### Data

For both training and external validation (testing) sets, patients from French centers were excluded for legal reasons (N = 118, not included in the numbers above). The training set comprised the FIR (NCT01846416)²⁵, POPLAR (NCT01903993)³ and BIRCH (NCT02031458)²⁶ phase 2 clinical trials. The test set was the atezolizumab arm of the OAK phase 3 trial (NCT02008227)²⁷ for individual predictions and additionally the docetaxel arm for trial predictions (Supplementary Figure 1). These studies were conducted in accordance with the Declaration of Helsinki after approval by institutional review boards or independent ethics committees. All patients provided written informed consent.

The outcome considered was overall survival (OS), defined as the time between treatment start and death or last follow-up, in which case the data was right-censored. The median follow-up was 35.2 months (95%CI:34.5–35.7) in the training set and 26.8 months (95%CI:26.3–27.5) in the test set.

### Preprocessing
## Baseline data

The baseline data consisted of 63 variables spanning demographic and biological data, clinical information and disease status (see Supplementary Figure 2–4 for a description of the main variables). PD-L1 expression on tumor cells was measured by immunohistochemistry or quantitative polymerase chain reaction, with four possible levels (0: < 1%; 1: ≥ 1%; 2: ≥ 5% and 3: ≥ 50%)³. We refer to the above-mentioned identifiers and references for further details on the other variables. Data were measured in accordance to the studies principles.

## Tumor and blood markers kinetics (TK and BK)

Patients with only one baseline SLD measurement and no SLD measurement during the treatment period were excluded (N = 110). For BK, first time points prior to treatment start were discarded. Then, four exclusion rules were established to identify anomalous data points: 1) values outside physiologically possible bounds, 2) duplicates, 3) values that abruptly went to an extreme out-of-range value between two measurements, 4) only the BK value at the closest time point to treatment initiation was kept. Eventually, in order to have sufficient data for Bayesian estimation with early data, patients with less than three observations before cycle 5 were removed. We refer to the supplementary methods for details.

## Nonlinear mixed-effects modeling

### Population approach

Statistical hierarchical nonlinear mixed-effects modeling (NLME) was used to implement a population approach²⁸ for the kinetic data and parameter estimation was conducted using the Monolix software²⁹. Mathematical details are given in the supplementary methods.

### Structural models

Following previous work, the TK structural model was assumed to be the sum of two exponentials¹⁹,³⁰:

$$y_j^i = \begin{cases} y_0^i e^{KG_i t} & t \le 0 \\ y_0^i (e^{-KS_i t} + e^{KG_i t} - 1) & t > 0 \end{cases}$$

where t = 0 corresponds to treatment initiation and y₀, KG and KS are three parameters, representing respectively the baseline value, growth and shrinkage rates. This model was also considered for BK, together with three other models: constant (yⱼ = αᵢ, ∀j), linear (yⱼ = αᵢ + βᵢtⱼ, ∀j) and hyperbolic (yⱼ = pᵢ + (eᵢt(qᵢ - pᵢ))³¹. Quantitative comparison of goodness-of-fit between models was assessed using the corrected Bayesian information criterion³².

### Identification of individual model-based parameters

The population parameters identified on the training set were used to define prior distributions of the TK and BK model parameters. These “training” priors were used for Bayesian estimation (maximum a posteriori estimate) of the individual TK and BK model parameters, not only for the training set but also for the test sets, in order to avoid leakage. To focus on the pure kinetic parameters, the model-estimated baseline parameters were not kept. We additionally considered the ratio of the model-predicted value at cycle 3 day 1 to the model-estimated baseline parameter. Altogether, there were three individual parameters for each marker: X KG, X KS and X ratio for X = TK, CRP, LDH and neutrophils; and albuminₚ, albuminᵢ and albumin ratio for albumin.

### Truncated data: individual-level

Individual-level truncated datasets were derived from the longitudinal TK and BK data by keeping data only up to: cycle 3 day 1 (C3D1, 1.5 months), C5D1 (3 months) and C10D1 (6.75 months). New training priors were estimated from each CXD1 training set. The resulting TK and BK truncated model parameter Y for marker X at cycle i were denoted by X Y, i (e.g., ldh KG, 5).
### Truncated data: study-level for trial prediction

128

Study-level truncated datasets were defined at the following on-study landmark times $lt$ after study initiation (first patient recruited): $lt = 10, 25$ and $60$ weeks.

129

130

Only the patients enrolled before this time and their data collected up to $lt$ was used. Note that here $t = 0$ corresponds to study initiation and thus patients in these datasets have varying follow-up duration (from 0 to $lt$), in contrast to individual-level truncated datasets.

131

132

133

### Machine learning

134

#### Data preparation

135

Missing values (1.6% total, maximum 12% in one variable) were imputed with the median for numeric variables and mode for categorical variables, learned on the training set, even when applied to the test set. All numeric variables were centered and scaled. Means and standard deviations were learned on the train and carried to the test set.

136

137

138

139

#### Models

140

Model elaboration and development was performed exclusively on the training set, using 10 folds cross-validation for predictive performances evaluation. Due to censoring in the data, survival models were used: proportional hazards Cox regression³³, extreme gradient boosting (XGB) with either Cox or accelerated failure time (AFT) models³⁴ and random survival forests (RSF)²³. Nested cross-validation with inner bagging in each 10-fold cross-validation outer loop was used to evaluate the benefit of tuning the hyperparameters³⁵. Improvement of the performances was negligible with hyperparameter tuning (Supplementary Figure 5). Therefore, we used the default values of the hyperparameters. For the final RSF model: number of trees *ntree* = 500, number of variables to possibly split at each node *mtry* = 5, minimum size of terminal node *nodesize* = 15, number of random splits for splitting a variable *nsplit* = 10.

141

142

143

144

145

146

147

148

149

150

#### Evaluation

151

Predictive performances were assessed for either discrimination (c-index and classification metrics at horizon times $τ$), calibration (calibration curves) or stratification (dichotomized KM survival curves). For each individual, the RSF model gives two prediction outputs: a scalar value termed "mortality" that we will refer to as "ML score", and time-dependent predicted survival curves²³. The former was used to compute the c-index using the *rcorr.cens* function of the *hmisc* R package³⁶,³⁷. For prediction of survival at a horizon time $τ$, we used the latter to compute model-predicted probabilities of death at $τ$. Unless otherwise specified, $τ = 12$ months. Survival-adapted metrics of predictive performance were used for sensitivity, specificity, area under the receiver-operator curve (ROC AUC) and negative and positive predictive value (NPV and PPV) to account for censoring³⁸,³⁹. For computation of accuracy, censored patients before $τ$ were discarded ($N = 17/396$ in the test set at 12 months). The optimal cut-points used for individual OS predictions on the test set were defined as the Kaplan-Meier estimated survival probability in the training set at $τ$ (0.257 at 6 months, 0.437 at 12 months, 0.634 at 24 months).

152

153

154

155

156

157

158

159

160

161

162

163

164

For patient stratification (dichotomized KM curves), the ML score was used, with models trained on the training set and predicted on the test set. In order to assess stratification abilities to capture the 20% of long-term survivors, cut-points were set at the 20th percentiles for each variable/score evaluated. This cut-point arbitrary definition was also motivated by the aim to ensure fair comparison between multiple parameters on the same data. Significance of differences in KM curves was established using the logrank test, and hazard ratios were computed using proportional hazards Cox regression.

165

166

167

168

169

170

171

### Variable selection and minimal signature

172

Variable selection was performed only for the BSL data. The method was based on two steps: 1) sorting the variables using least absolute shrinkage and selection operator (LASSO)⁴⁰ and 2) building RSF incremental models including increasing numbers of variables. LASSO sorting was defined

173

174

175
as taking the coefficients gradually becoming non-zero during likelihood maximization when the regularization parameter decreases. The minimal signature was defined as the minimal set of variables able to achieve a c-index larger than 0.75 and an AUC larger than 0.8, with the addition of 4 well-established prognosis.

## Survival simulations and computation of predicted HRs

For each patient *i*, one output of the kML model is a survival curve $S^i(t)$. This gives the cumulative distribution function $1 - S^i(t)$ of the random variable $T^i$ of the time to death for patient *i*, which was used to simulate 100 replicates of $T^i$. Pooling all patients together, we thus obtained 100 replicates of {$T^i, ATZ, T^j,DTX$} for *i* and *j* being the patient indices within the ATZ and docetaxel arms, respectively. Each replicate then led to 1) a predicted survival curve in each arm and 2) a Cox proportional hazard HR between the two arms. Taking the mean and the 5th and 95th percentiles over all replicates yielded the reported point estimate and corresponding 95% prediction interval. The same procedure was used for study-truncated data.

## Data Availability

Qualified researchers may request access to individual patient level data through the clinical study data request platform (https://vivli.org/). Further details on Roche's criteria for eligible studies are available here (https://vivli.org/members/ourmembers/). For further details on Roche's Global Policy on the Sharing of Clinical Information and how to request access to related clinical study documents, see here https://www.roche.com/innovation/process/clinical-trials/data-sharing/.

## Code availability

Algorithms used for data analysis are all publicly available from the indicated libraries and references in the Methods section.

# Results

## Data

The data consisted of advanced NSCLC patients enrolled in ATZ trials (N = 1936, Figure 1B and Supplementary Figure 1). Three ATZ phase 2 trials were pooled into a training dataset<sup>3,25,26</sup> (N = 862). The external validation (test) set comprised data from the ATZ arm (N = 553) of the OAK phase 3 trial<sup>27</sup>. For trial outcome prediction the docetaxel arm (N = 521) was added as an additional test set.

Variables comprised baseline (pre-treatment) and longitudinal (on-treatment) data (Figure 1A). The former included: patients and disease characteristics (*p* = 63 variables, 43 numeric and 20 categorical, denoted BSL) and transcriptomic ("RNAseq", *p* = 58,311 transcripts) data. The latter included: longitudinal investigator-assessed sum of largest diameters (SLD) of lesions as per the RECIST criteria<sup>41</sup>, denoted by tumor kinetics (TK, *k* = 5,473/3,015 time points in the train/test sets, respectively, median 5/4 data points per patient, range 2/2–24/20); and longitudinal measurements of four blood markers (albumin, C-reactive protein (CRP), lactate dehydrogenase (LDH) and neutrophils), denoted together as blood markers kinetics (BK, *k* = 60,779/38,460 data points, median 11–7–11–11/9–9–9–10 data points per patient, range 3–3–3–3/3–3–3–3–60–63–63–78/82–47–77–89 for albumin–CRP–LDH-neutrophils in the train/test sets, respectively).

## Nonlinear mixed-effects modeling (NLME) of longitudinal markers

We first developed NLME models for the longitudinal data (Figure 1B). The TK structural model was the sum of an increasing and a decreasing exponential function (double exponential model)<sup>30</sup>. It was able to accurately describe the training data with no goodness-of-fit misspecification (Figure 2A and Supplementary Figure 6). Population parameters were estimated with good accuracy (all relative standard errors smaller than 9%, Table 1).
<figure><img src="image_1.png" /><figcaption>**Figure 1. Study schematic A.** Baseline and longitudinal data were combined into a machine learning algorithm in order to predict individual survival prognosis. Longitudinal data were modelled using nonlinear mixed-effects modelling, whereas machine learning-based feature selection was applied to the baseline data to derive a minimal signature. Tumor kinetics and biological kinetics parameters were combined with the minimal signature to predict survival. Predictive performances were assessed using survival metrics (c-index and survival at horizon times).</figcaption></figure>

**B. Algorithm** used to develop the model on the train data and carry it to the test set for external validation. Each step — preprocess, learning of the Bayesian priors, dimensionality reduction, feature selection, choice, tuning and training of the machine learning algorithm — were calibrated on the training set and then applied to the test set.

TK: tumor kinetics; BK: blood markers kinetics; ML: machine learning; NLME: nonlinear mixed-effects modelling
<table><caption>**Table 1.** Parameters from nonlinear mixed-effects modeling of tumor and blood marker kinetics</caption><thead><tr><th></th><th>K</th><th>CRP</th><th>LDH</th><th>Neutrophils</th><th>Albumin</th><th></th></tr></thead><tbody><tr><td>KG<sub>pop</sub> (week<sup>-1</sup>)</td><td>0.004 92</td><td>(6.80)</td><td>0.008 14</td><td>( 9.38)</td><td>0.002 38</td><td>(10.48)</td><td>0.004 36</td><td>( 8.69)</td><td>p<sub>pop</sub> (g/l)</td><td>29.4</td><td>(3.82)</td></tr><tr><td>KG<sub>pop</sub> (week<sup>-1</sup>)</td><td>0.007 78</td><td>(8.22)</td><td>0.0137</td><td>(14.14)</td><td>0.001 84</td><td>(13.73)</td><td>0.000 987</td><td>(21.16)</td><td>l<sub>pop</sub> (log (day))</td><td>8.09</td><td>(2.74)</td></tr><tr><td>ω<sub>KG</sub></td><td>1.36</td><td>(3.80)</td><td>1.61</td><td>( 4.25)</td><td>1.55</td><td>( 5.36)</td><td>1.41</td><td>( 4.48)</td><td>ω<sub>p</sub></td><td>0.476</td><td>(7.48)</td></tr><tr><td>ω<sub>KS</sub></td><td>1.41</td><td>(4.66)</td><td>1.81</td><td>( 6.29)</td><td>1.92</td><td>( 5.34)</td><td>2.46</td><td>( 5.82)</td><td>ω<sub>l</sub></td><td>0.359</td><td>(6.42)</td></tr><tr><td>error<sup>1</sup></td><td>6.82</td><td>(1.15)</td><td>0.559</td><td>( 1.23)</td><td>0.138</td><td>( 0.79)</td><td>0.207</td><td>( 0.82)</td><td>error<sup>1</sup></td><td>0.0549</td><td>(0.77)</td></tr></tbody></table>

Parameter value (relative standard error (%)). TK: constant error, others: proportional error. CRP: C-reactive protein; LDH: lactate dehydrogenase.

To analyze the BK data, we first investigated whether significant kinetic patterns could be observed beyond random noise (due to, e.g., measurement errors, see raw data in Supplementary Figures 7–10). The latter was considered as the null hypothesis, described by a constant model. It was tested against three alternative empiric models: linear, hyperbolic (monotonous but non-linear and saturating) and double-exponential (nonlinear and non-monotonous). For all four BKs, we found significant kinetics compared with the constant model, as shown by lower corrected Bayesian information criterion and relative error between model fits and data (Supplementary Figure 11). The best descriptive models were hyperbolic for albumin and double-exponential for the other BKs. Individual fits to patient kinetics with the best models showed substantial descriptive power (Figure 2A), which was confirmed by data versus model fits plots (Supplementary Figures 12–15). Parametric identifiability of population parameters was excellent for all models (Table 1).

We further assessed the stratification value of the individual model-based kinetic marker for OS prognosis (Figure 2B). The TK parameter KG (growth rate) exhibited good stratifying ability (HR = 4.39 (2.8–6.89)), which was similar to the $CRP_{KG}$ parameter (HR = 4.37 (2.76–6.91)). Ranked by HR importance; (controlled by the 20th percentile definition of the cut-point, see methods), the following four best parameters were $albumin_p$ (HR = 3.17 (2.11–4.78)), neutrophils$_{KG}$ (HR = 3.07 (2.04–4.63)), neutrophils$_{KS}$ (HR = 2.33 (1.6–3.39)) and TK$_{KS}$ (HR = 2.02 (1.42–2.89)). All kinetic parameters carried substantial prognostic power ($p < 0.0001$, log rank test).

For TK and BKs we complemented the initial model parameters with an additional metric that was considered valuable for early prediction: the model-predicted ratio of change over baseline at cycle 3 day 1.

## Survival prediction using kinetics-machine learning (kML): model development

Four feature sets resulted from the analysis above: BSL, RNAseq, TK and BK (Figure 1A). The development of a kinetics-machine learning (kML) comprised two main steps: choice of the algorithm and derivation of a minimal signature (Figure 1B). The first was achieved by benchmarking four models that used all variables ($p = 119$, $N = 553$). The random survival forest (RSF) model found to exhibit the best performances (Supplementary Figure 5) and was thus selected. Notably, we found significantly better predictive performances of RSF over a classical Cox proportional hazard regression model ($p = 0.0006$).

Feature selection on BSL variables was performed building incremental RSF models based on LASSO importance-sorted variables (Figure 3A). The model using all of them achieved the best score. Nevertheless, keeping in mind the objective to ultimately support decision making and patient stratification, a minimal (11 features), near-optimal, set of BSL variables was selected and denoted mBSL. It was defined as the first seven variables reaching the plateau (CRP, heart rate, neutrophils to lymphocytes ratio, neutrophils, lymphocytes to leukocytes ratio, liver metastases and ECOG score), complemented with four variables with established prognostic or predictive value and available in routine care: PD-L1 expression (50% cut-off)³, hemoglobin⁴², SLD²² and LDH⁴³,⁴⁴.

Applying stringent criteria to the RNAseq data (see supplementary methods), we selected 167 transcripts as candidates for final variable selection using Bolasso regression model to identify the
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Goodness-of-fit metrics and plots of dynamic BK models A. Representative individual fits for the TK and BK best empirical models showing non-trivial kinetic parameters well captured by the dynamic models. Survival is indicated by a vertical line (solid = death, dashed = censored). B. Stratified Kaplan-Meier curves at the 20th percentile level on the test set, for TK and BK model-based parameters. Missing values were removed in this univariable analysis, explaining the difference of initial number of patients for albumin that had 9 patients in this case. CR: complete response; PR: partial response; SD: stable disease; PD: progressive disease.</figcaption></figure>
<table><caption>**Table 2.** Contingency table for OS prediction at 12 months</caption><thead><tr><th rowspan="2">MODEL</th><th rowspan="2"></th><th colspan="3">TRUTH</th><th rowspan="2">Total</th></tr><tr><th>Alive (0)</th><th>dead (1)</th><th>Total</th></tr></thead><tbody><tr><td></td><td>Alive (-)</td><td>182</td><td>30</td><td>212</td><td>(58.7%)</td></tr><tr><td></td><td>Dead (+)</td><td>48</td><td>101</td><td>149</td><td>(41.3%)</td></tr><tr><td></td><td>Total</td><td>230</td><td>(63.7%)</td><td>131</td><td>(36.3%)</td></tr></tbody></table>

Note: 16/377 censored patients with survival time ≤ 12 months removed for computation of accuracy. sensitivity, specificity, PPV and NPV don't correspond exactly to the numbers because they are computed from KM estimate, thus adjusting for censoring bias.

optimal set of predictors⁴⁵. Finally, we ended up with 52 RNAseq variables that corresponded to the highest average c-index of 0.64.

We then compared the cross-validated c-index of each feature set on the train data (Figure 3B). Because of negligible discrimination performances (c-index = 0.62 ± 0.050) and non-systematic availability of those data, the RNAseq set was removed from the model. The selected set of clinical data at baseline (mBSL) exhibited moderate discrimination performances (c-index = 0.710 ± 0.038), which was slightly outperformed by the TK set (c-index = 0.723 ± 0.025). Interestingly, the BK set significantly outperformed both baseline clinical and TK (c-index = 0.793±0.038, *p* = 0.0004 and 0.0005 respectively, Student's t-test). Jointly, mBSL, TK and BK performed significantly better than any feature set alone (c-index = 0.824 ± 0.050, *p* = 0.00007, 0.0002 and 0.055), as well as any combination of two sets among the three (mBSL + TK: c-index = 0.77 ± 0.026, mBSL + BK: c-index = 0.81 ± 0.027, TK + BK: c-index = 0.80 ± 0.049). The resulting model combining mBSL, TK and BK was denoted kML (kinetics-machine learning).

During cross-validation on the training set, kML exhibited excellent predictive performances across multiple metrics, with minimal between-folds variability (AUC = 0.919 ± 0.056, accuracy = 0.873 ± 0.052, Figure 3C).

## External validation

The predictive performance of the final kML model (mBSL, TK and BK) was assessed on the ATZ test set (377 patients). At the population level, the model-predicted survival curve was in excellent agreement with the observed data (Figure 4A). Notably, the prediction interval from the model was narrow, indicating high precision. At the individual level, consistent with the cross-validation results, substantial discrimination performances were observed (c-index = 0.790, accuracy and AUC for 12-months survival probability 0.787 and 0.874, respectively, Figure 4B). All classification metrics for prediction of survival at 12 months were high (≥ 0.78), except PPV, indicating worse ability to predict death than survival. Although smaller, they were similar to the cross-validation results.

In addition, calibration curves revealed good performance, at multiple horizon times (Figure 4C). Model-predicted probabilities were concordant with the observed KM estimates of the survival probabilities, over the entire range of the binned predicted probabilities. This is further illustrated by the contingency Table 2. For instance, among 212 patients predicted to be alive at 12 months, 182 (85.8%) were actually alive. Predictive AUC was good at other horizon times (0.846 and 0.910 at 6 and 24 months, respectively, Supplementary Figure 16). However, PPV and sensitivity were very low at 6 months.

Notably, the kML mortality score derived from the model and learned on the training set was able to accurately stratify OS in the test set (HR = 25.2 (10.4–61.3), *p* < 0.0001, Figure 4D), indicating excellent ability to identify the 20% of long-term survivors. It outperformed all single kinetic markers (Figure 2C).

Variables importance was assessed by running a post-hoc multivariable Cox regression (Figure
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Minimal baseline (mBSL) signature and kinetics-ML (kML) model A. Cross-validated (CV) performance scores on the training set (c-index and AUC, mean ± standard deviation) for incremental random survival forest (RSF) models using an increasing number of baseline clinical and biological variables sorted by LASSO importance. The dashed blue line shows the minimal number of variables reaching the plateau. Blue-colored variables correspond to the minimal clinical signature (mBSL). B. Comparative CV c-indices of RSF models based either on RNAseq, mBSL, TK, BK and mBSL + TK + BK (final model, kML) variables showing increased predictive performances over baseline when using model-based parameters of kinetic markers. Numbers on the bars indicate the number of variables. C. CV performances of the kML model for discrimination (c-index) and classification (survival prediction at 12-months OS).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4: Predictive performances of kML on the ATZ test set**</figcaption></figure>

**Figure 4.** Predictive performances of kML on the ATZ test set **A.** Comparison of the population-level survival curves between the data (KM estimator) and the model prediction. **B.** Scores of discrimination metrics. Classification metrics were computed for prediction of OS at 12 months. **C.** Calibration curves at 6, 12 and 24 months, showing the observed survival probabilities (with KM 95% confidence interval) versus the predicted ones in 10 bins corresponding to the model-predicted survival probability deciles. Dashed line is the identity. **D.** Dichotomized KM survival curves based on the ML model-predicted score (high versus low), at the 20th percentile cut-off. **E.** Variables importance (multivariable hazard ratios) in the full time-course kML model.

4F). Interestingly, the top two variables were BKs (CRP~KG~ and CRP ratio C3). In addition, TK and BK made up for six out of the seven top important features and were found more important than PD-L1.

Given the large sample size of our data, we further assessed the model performances when trained on smaller data sets (Supplementary Figure 17). The learning curve revealed that approximately 200 patients were necessary to reach similar performance to the ones obtained with the full training set (N = 533), for both cross-validation and external validation on the test set (c-index = 0.82 ± 0.056 vs c-index = 0.82 ± 0.050 in cross-validation, 0.78 vs 0.79 on the test set, models trained with 200 vs 533 patients, respectively). Trained with only 60 patients, kML reached already good performances (c-index = 0.76 ± 0.15 and 0.74 in cross-validation and test, respectively).

Together, these results demonstrate important predictive performances of overall survival following ATZ treatment using the kML model.

**Application to individual survival prognosis from early on-treatment data**

Results above required full on-treatment time-course data to compute TK and BK markers, thus cannot be used to make early predictions. To investigate the operational applicability of our methodology, data from the test set were truncated at the beginning of treatment cycles number 3, 5 and 10, respectively corresponding to 1.5, 3 and 6.75 months. We found that integrating longer on-treatment data in kML, the predictive performances steadily increased (Figure 5A and Supplementary Figure 18). Using the baseline variables only (mBSL), the stratification ability was significant but moderate (HR = 1.74 (1.24–2.46), *p* = 0.0014, Figure 5B). In contrast, kML exhibited increasing stratification ability from data at 1.5 months (HR = 2.19 (1.53–3.12), *p* < 0.0001), 3 months (HR = 3.51
<figure><img src="image_5.png" /><figcaption>**Figure 5. Predictive value of kML from cycle-truncated data A.** Predictive power (c-index) of ML models using baseline (BSL) or truncated data at 1.5, 3 and 6.8 months as well as the full time-course. **B.** Stratified KM survival curves using a RSF model trained on the minimal baseline (mBSL) variables. **C.** Stratified KM survival curves using kML from 1.5 months (2 cycles), 3 months (4 cycles) and 6.8 months (9 cycles) truncated data. Truncation time is indicated by the vertical line.</figcaption></figure>

TK: tumor kinetics; BK: biological kinetics; LDH: lactate dehydrogenase; CRP: C-reactive protein.

(2.33 – 5.3), *p* < 0.0001) and 6.8 months (HR = 5.01 (3.16 – 7.95), *p* < 0.0001), see Figure 5C.

Further investigation of the predictive performances of individual kinetic markers revealed that TK parameters were the most informative at 6 weeks (1.5 months, first imaging assessment). Adding BKs to TKs brought additional predictive value starting at 3 months, and BKs outperformed TK from 6.75 months on (Supplementary Figure 19A). Among BKs, neutrophils kinetics appeared to be the most predictive, followed by CRP, albumin and LDH. However, the combined BK signature outperformed each individual BK, indicating that their collective predictive capabilities were not driven by any single biomarker alone.

Interestingly, the most important variable at 1.5 months was a kinetic one, TK ratio C3 with following variables being from mBSL (e.g., liver metastases, PDL1 and ECOG). When more on-treatment variables become available, this shifted to TK and BK (TK ratio C3, $TK_{KS}$, $TK_{KG}$, $CRP_{KG}$, $LDH_{KG}$), see Supplementary Figure 19B.

Application to clinical trial outcome prediction from early on-study data

The kML model can also be applied for the prediction of the outcome of a clinical trial (survival curves and associated hazard ratio), from early on-study data. We performed on-study runca-tions on the test set based on a number of weeks after the date of the first patient recruited (see methods). Here, we applied the model to predict not only patients receiving ATZ, but also doc-etaxel (Figure 1B). Predictions of the kML model applied to each arm yielded very accurate results when using data from the entire study (predicted HR = 0.784 (0.7 – 0.842)), versus data HR = 0.778 (0.65 – 0.931), Figure 6A–B). Notably, the model prediction intervals were narrower than the data Kaplan-Meier confidence intervals, probably because the kML-trained model incorporates the in-formation from the three phase 2 trials. Using only early data, the model was already able to detect
<figure><img src="image_6.png" /><figcaption>Figure 6: Use of kML for early prediction of the outcome of a clinical trial</figcaption></figure>

**Figure 6.** Use of kML for early-prediction of the outcome of a clinical trial A. Survival curves model-based predictions and prediction intervals versus actual data from on-study data at multiple horizon times after study initiation. Note that the model is able to predict full survival curves even if based on early kinetics. B. Compared data and kML-predicted hazard ratios. C. Description of hazard ratios, number of patients and number of data points available in each arm, at the landmark on-study time points. PI: prediction interval, CI: confidence interval, DTX: docetaxel arm, ATZ: atezolizumab arm.

a (non-significant) tendency at 10 weeks, with only 23 and 30 patients in each arm, and very short follow-up. Starting from data available at 25 weeks (6.25 months), the model correctly predicted a positive outcome of the study, with a 95% prediction interval of the HR below 1. Of note, the available data at this time (dashed lines, Figure 6A and red HR CIs in Figure 6B) was far from being conclusive. The model prediction was stable from 25 weeks on whereas the OS data only exhibited significant HR starting from 60 weeks and required more than 300 patients in each arm to be conclusive.

Discussion

Blood markers from hematology and biochemistry are routinely collected during clinical care or drug trials. They are cost-effective and easily obtained both before and during treatment. There is limited exploration regarding the predictive capabilities of the kinetics of such data. Combining BSL variables with on-treatment data (TK and BK), we addressed this question using a novel hybrid NLME-ML methodology. The resulted kML model demonstrated excellent predictive perfor-
mances for OS in two aspects: 1) patient-level predictions (discrimination, calibration and patient stratification) and 2) trial-level predictions. The kML model outperformed current state-of-the-art methods based on either baseline or on-treatment data alone, utilizing only routine clinical information, with a c-index of 0.79 and an accuracy of 78% for prediction of 12-month survival, on the test dataset. Overall, kML incorporates 26 features, out of which 15 features require monitoring five quantities over time (tumor size, albumin, CRP, LDH and neutrophils).

Regarding baseline markers, the predictive value of PD-L1 expression, commonly used in clinical care, is controversial<sup>9,10</sup>. Previous studies reported an AUC for durable response of 0.601 and a PFS HR of 1.90 (PD-L1 ≥ 1% vs 0%)<sup>8</sup>. Baseline tumor mutational burden showed similar predictive value initially (AUC = 0.646)<sup>11</sup>, but led to disappointing results in a recent prospective study<sup>46</sup>. Baseline blood counts were previously reported to predict overall survival<sup>43,47-49</sup> and treatment response (AUC = 0.74)<sup>42</sup>. The ROPRO score, derived from a large pan-cancer cohort and incorporating baseline clinical and biological data (27 variables) achieved a c-index of 0.69 and a 3-months AUC of 0.743 for prediction of survival in the OAK clinical trial<sup>50</sup>. Here, we confirmed these findings and established a minimal signature of such data composed of only 11 variables (CRP, heart rate, neutrophils to lymphocytes ratio, neutrophils, lymphocytes to leukocytes ratio, liver metastases, ECOG, PD-L1 ≥ 50%, hemoglobin, SLD and LDH), yet with similar predictive performances (c-index = 0.678) and significant stratification ability (HR = 1.74, p = 0.0014). Altogether, our kML model demonstrated substantially better predictive performances than these baseline models.

We further confirmed the established predictive value of TK model-based parameters<sup>19,20</sup>. Blood- or serum-derived longitudinal markers kinetics have to date rarely been modeled. Gavrilov et al. proposed to model NLR kinetics and demonstrated improved OS predictions over TK alone<sup>31</sup>. Here we extended to four BKs: albumin, CRP, LDH and neutrophils. This choice was not only motivated by observed statistical associations, but also from biological considerations. Albumin is associated with nutritional status (cachexic state) and is known to evolve with time in responders. CRP is a marker of systemic inflammation<sup>44</sup>. Increased CRP, decreased albumin level, and increased CRP/albumin ratio have been reported to be associated with poor survival<sup>51</sup>. Neutrophils play a role in inflammation by promoting a favorable microenvironment for cancer cell growth and spread, and activation of carcinogenic signaling pathways<sup>52</sup>. Elevated LDH levels are a marker of cancer cells turnover rate, and LDH has a potential role for prediction of potential invisible metastases<sup>44</sup>. We found that all these markers had non-trivial on-treatment kinetics. However, data fits were not perfect, possibly due to the simplicity and empiric nature of the models we used. Further mechanistic modeling of the joint kinetics of BKs and TK could bring relevant biological information and yield more accurate predictive parameters. We found that all four BKs were contributive to the model and that, combined, they outperformed TK performances.

We analyzed the RNAseq data using standard methods and found only negligible predictive performances. Such result could be explained by the fact that the tissue of origin that was used was heterogeneous across the patients (primary tumor or metastasis), was limited to a local area of the tumor, and could come from tissue sampled long before treatment initiation. Given that our main objective was to derive a predictive model from markers available in routine practice, we excluded it from our minimal signature. A refined analysis, especially focusing on immune-based signatures, could improve our results<sup>5</sup>.

Machine learning models, although increasingly used in pharmacological studies—including recently for TK-OS modeling and variable selection<sup>22,53</sup>—have yet rarely been rigorously compared to classical statistical models<sup>24</sup>. Here, such comparison revealed significantly better performance of the nonlinear random survival forest RSF model compared to the linear proportional hazards Cox model. In our approach, we did not use the propagation of standard statistical quantification of the parameters’ estimates uncertainty to evaluate the accuracy of the model predictions. Rather, we relied on the RSF-outputted individual survival curves to sample virtual individuals and compute prediction intervals.

A drawback of classical TK-OS studies is that they make use of the full observed kinetics to
predict overall survival, which can lead to time-dependent covariate bias<sup>54</sup> and limit their practical applicability at bedside. We used individual-truncated data sets and found that kML was already improving predictions over mBSL using data at 1.5 months, which corresponds to the first imaging assessment of the treatment effect. At later times, stratification abilities increased to highly significant levels (e.g., HR = 5 at 6.8 months).

A strength of our study is that we relied on well-curated data with high number of patients from clinical trials. However, when extrapolating to other settings — earlier trial phases, real-world data — limited number of patients might be available. Yet, we found that using only 60 patients to train kML was sufficient to reach near-optimal performances.

Not only kML has value for personalized health care, but it also revealed useful for prediction of a phase 3 trial using early on-study data. Our model was able to predict the study's positive outcome with data at ~ 6 months, versus 10 months using TK only<sup>19</sup>. Relying on the data alone, such positive outcome was only detectable at 15 months. These results could have important implications for drug development as they could inform earlier on go/no-go decisions. Consequently, this could allow to detect futility more easily and more rapidly during clinical trials, allowing to avoid treating patients with an inefficient investigational treatment and to reassign funds and energy to other researches. Of note, in a recent evaluation based on resampling the first-line NSCLC ATZ study IMpower150 to mimic small, short follow up early Phase Ib studies, TK model-based metrics had better operating characteristics to predict Phase III success compared with RECIST endpoints ORR and PFS<sup>55</sup>. Extension of such results with the addition of BKs is thus a promising line of research. In addition, kML, trained on ATZ data, yielded excellent predictive abilities for the docetaxel (control) arm. This suggests that the relationships between TK / BK and OS might be drug-independent. In turn, this opens future perspectives in terms of testing kML on drugs with different mechanism of action, or combinations.

Further avenues of research comprise the development of integrative models from advanced multi-modal data such as the one collected during the PIONeeR clinical study (NCT03833440)⁵⁶,⁵⁷ that include quantitative image analysis from multiplex immune-histochemistry, genomic and transcriptomic data, biological and clinical markers. In addition, mechanistic modeling of quantitative and physiologically meaningful longitudinal data (immune-monitoring, vasculo-monitoring, circulating DNA¹²,⁵⁸⁻⁶⁰, soluble factors⁶¹, pharmacokinetics, TK and a large number of BKs from either hematology or biochemistry) paves the way to an improved understanding and prediction of mechanisms of relapse to ICI⁶². Furthermore, the predictive abilities of kML — at the both individual and study levels — should be evaluated in model-based prospective trials⁶³.

In conclusion, our study shows that integrating model-based on-treatment dynamic data from routine biological markers shows great promise for both personalized health care and early prediction of the outcome of clinical trials during drug development.

## References

1. Bray, F., Ferlay, J., et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. *CA: A Cancer Journal for Clinicians*, 394–424. ISSN: 1542-4863. doi:10.3322/caac.21492 (2018).

2. Duma, N., Santana-Davila, R. & Molina, J. R. Non-Small Cell Lung Cancer: Epidemiology, Screening, Diagnosis, and Treatment. *Mayo Clinic Proceedings*, 1623–1640. ISSN: 0025-6196, 1942-5546. doi:10.1016/j.mayocp.2019.01.013 (2019).

3. Fehrenbacher, L., Spira, A., et al. Atezolizumab versus docetaxel for patients with previously treated non-small-cell lung cancer (POPLAR): a multicentre, open-label, phase 2 randomised controlled trial. *The Lancet*, 1837–1846. ISSN: 0140-6736, 1474-547X. doi:10.1016/S0140-6736(16)00587-0 (2016).
451

452

453

4. Grant, M. J., Herbst, R. S. & Goldberg, S. B. Selecting the optimal immunotherapy regimen in driver-negative metastatic NSCLC. *Nature Reviews Clinical Oncology*, 625-644. ISSN: 1759-4782. doi:10.1038/s41571-021-00520-1 (2021).

454

455

456

5. Camidge, D. R., Doebele, R. C. & Kerr, K. M. Comparing and contrasting predictive biomarkers for immunotherapy and targeted therapy of NSCLC. *Nature Reviews Clinical Oncology*, 341-355. ISSN: 1759-4782. doi:10.1038/s41571-019-0173-9 (2019).

457

458

6. Hutchinson, L. & Kirk, R. High drug attrition rates—where are we going wrong? *Nature Reviews Clinical Oncology*, 189-190. ISSN: 1759-4782. doi:10.1038/nrclinonc.2011.34 (2011).

459

460

461

462

7. Hua, T., Gao, Y., Zhang, R., Wei, Y. & Chen, F. Validating ORR and PFS as surrogate endpoints in phase II and III clinical trials for NSCLC patients: difference exists in the strength of surrogacy in various trial settings. *BMC Cancer*, 1022. ISSN: 1471-2407. doi:10.1186/s12885-022-10046-z (2022).

463

464

465

466

8. Rizvi, H., Sanchez-Vega, F., et al. Molecular Determinants of Response to Anti-Programmed Cell Death (PD)-1 and Anti-Programmed Death-Ligand 1 (PD-L1) Blockade in Patients With Non-Small-Cell Lung Cancer Profiled With Targeted Next-Generation Sequencing. *Journal of Clinical Oncology*. doi:10.1200/JCO.2017.75.3384 (2018).

467

468

469

9. Doroshow, D. B., Bhalla, S., et al. PD-L1 as a biomarker of response to immune-checkpoint inhibitors. *Nature Reviews Clinical Oncology*, 345-362. ISSN: 1759-4782. doi:10.1038/s41571-021-00473-5 (2021).

470

471

472

10. So, W. V., Dejardin, D., Rossmann, E. & Charo, J. Predictive biomarkers for PD-1/PD-L1 checkpoint inhibitor response in NSCLC: an analysis of clinical trial and real-world data. *Journal for Immunotherapy of Cancer*, e006464. ISSN: 2051-1426. doi:10.1136/jitc-2022-006464 (2023).

473

474

475

11. Hellmann, M. D., Ciuleanu, T.-E., et al. Nivolumab plus Ipilimumab in Lung Cancer with a High Tumor Mutational Burden. *New England Journal of Medicine*, 2093-2104. doi:10.1056/NEJMoa1801946 (2018).

476

477

478

12. Gandara, D. R., Paul, S. M., et al. Blood-based tumor mutational burden as a predictor of clinical benefit in non-small-cell lung cancer patients treated with atezolizumab. *Nature Medicine*, 1441-1448. ISSN: 1546-170X. doi:10.1038/s41591-018-0134-3 (2018).

479

480

13. Cristescu, R., Mogg, R., et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. *Science*, eaar3593. doi:10.1126/science.aar3593 (2018).

481

482

14. Sankar, K., Ye, J. C., et al. The role of biomarkers in personalized immunotherapy. *Biomarker Research*, 32. ISSN: 2050-7771. doi:10.1186/s40364-022-00378-0 (2022).

483

484

15. Acosta, J. N., Falcone, G. J., Rajpurkar, P. & Topol, E. J. Multimodal biomedical AI. *Nature Medicine*, 1773-1784. ISSN: 1546-170X. doi:10.1038/s41591-022-01981-2 (2022).

485

486

487

16. Kurtz, D. M., Esfahani, M. S., et al. Dynamic Risk Profiling Using Serial Tumor Biomarkers for Personalized Outcome Prediction. *Cell*, 699-713.e19. ISSN: 1097-4172. doi:10.1016/j.cell.2019.06.011 (2019).

488

489

17. Bonate, P. L. *Pharmacokinetic-Pharmacodynamic Modeling and Simulation* 2nd ed. 2011. ISBN: 978-1-4419-9484-4 (Springer-Verlag New York Inc., New York, 2011).

490

491

492

18. Claret, L., Girard, P., et al. Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics. *J Clin Oncol*, 4103-4108. doi:10.1200/JCO.2008.21.0807 (2009).

493

494

495

19. Claret, L., Jin, J. Y., 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. *Clin Cancer Res*, 3292-3298. doi:10.1158/1078-0432.CCR-17-3662 (2018).
20. Chan, P., Marchand, M., et al. Prediction of overall survival in patients across solid tumors following atezolizumab treatments: A tumor growth inhibition-overall survival modeling framework. *CPT: pharmacometrics & systems pharmacology*, 1171-1182. ISSN: 2163-8306. doi:10.1002/psp4.12686 (2021).

21. Benzekry, S. Artificial intelligence and mechanistic modeling for clinical decision making in oncology. *Clinical Pharmacology & Therapeutics*, 471-486. ISSN: 1532-6535. doi:10.1002/cpt.1951 (2020).

22. Chan, P., Zhou, X., et al. Application of Machine Learning for Tumor Growth Inhibition - Overall Survival Modeling Platform. *CPT: pharmacometrics & systems pharmacology*, 59-66. ISSN: 2163-8306. doi:10.1002/psp4.12576 (2021).

23. Ishwaran, H., Kogalur, U. B., Blackstone, E. H. & Lauer, M. S. Random survival forests. *The Annals of Applied Statistics*, 841-860. ISSN: 1932-6157, 1941-7330. doi:10.1214/08-AOAS169 (2008).

24. Christodoulou, E., Ma, J., et al. A systematic review shows no performance benefit of machine learning over logistic regression for clinical prediction models. *Journal of Clinical Epidemiology*, 12-22. ISSN: 0895-4356. doi:10.1016/j.jclinepi.2019.02.004 (2019).

25. Spigel, D. R., Chaft, J. E., et al. FIR: Efficacy, Safety, and Biomarker Analysis of a Phase II Open-Label Study of Atezolizumab in PD-L1-Selected Patients With NSCLC. *Journal of Thoracic Oncology*, 1733-1742. ISSN: 1556-0864. doi:10.1016/j.jtho.2018.05.004 (2018).

26. Peters, S., Gettinger, S., et al. Phase II Trial of Atezolizumab As First-Line or Subsequent Therapy for Patients With Programmed Death-Ligand 1-Selected Advanced Non-Small-Cell Lung Cancer (BIRCH). *Journal of Clinical Oncology: Official Journal of the American Society of Clinical Oncology*, 2781-2789. ISSN: 1527-7755. doi:10.1200/JCO.2016.71.9476 (2017).

27. Rittmeyer, A., Barlesi, F., et al. Atezolizumab versus docetaxel in patients with previously treated non-small-cell lung cancer (OAK): a phase 3, open-label, multicentre randomised controlled trial. *The Lancet*, 255-265. ISSN: 0140-6736, 1474-547X. doi:10.1016/S0140-6736(16)32517-X (2017).

28. Lavielle, M. *Mixed Effects Models for the Population Approach* ISBN: 1-4822-2650-2 (CRC Press, 2014).

29. Lixoft. *Monolix version 2020R1*. Antony, France, 2020.

30. Stein, W. D., Figg, 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*, 1046-1054. doi:10.1634/theoncologist.2008-0075 (2008).

31. Gavrilov, S., Zhudenkov, K., et al. Longitudinal Tumor Size and Neutrophil-to-Lymphocyte Ratio Are Prognostic Biomarkers for Overall Survival in Patients With Advanced Non-Small Cell Lung Cancer Treated With Durvalumab. *CPT: pharmacometrics & systems pharmacology*, 67-74. ISSN: 2163-8306. doi:10.1002/psp4.12578 (2021).

32. Delattre, M., Lavielle, M. & Poursat, M.-A. A note on BIC in mixed-effects models. *Electronic Journal of Statistics*, 456-475. ISSN: 1935-7524, 1935-7524. doi:10.1214/14-EJS890 (2014).

33. Cox, D. R. Regression Models and Life-Tables. *Journal of the Royal Statistical Society. Series B (Methodological)*, 187-220. ISSN: 0035-9246 (1972).

34. Chen, T. & Guestrin, C. XGBoost: A Scalable Tree Boosting System in *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining* (Association for Computing Machinery, San Francisco, California, USA, 2016), 785-794. ISBN: 9781450342322. doi:10.1145/2939672.2939785.

35. Cawley, G. C. & Talbot, N. L. C. On over-fitting in model selection and subsequent selection bias in performance evaluation. *Journal of Machine Learning Research*, 2079-2107 (2010).
543

544

36. Harrell, F. E., Califf, R. M., Pryor, D. B., Lee, K. L. & Rosati, R. A. Evaluating the yield of medical tests. *JAMA*, 2543-2546. ISSN: 0098-7484 (1982).

545

546

37. Harrell, F. E. Hmisc: Harrell miscellaneous. *R package*. https://CRAN.R-project.org/package=Hmisc (2022).

547

548

549

38. Heagerty, P. J., Lumley, T. & Pepe, M. S. Time-dependent ROC curves for censored survival data and a diagnostic marker. *Biometrics*, 337-344. ISSN: 0006-341X. doi:10.1111/j.0006-341X.2000.00337.x (2000).

550

551

39. Heagerty, P. J. & Saha-Chaudhuri, p. b. P. survivalROC: Time-dependent ROC curve estimation from censored survival data tech. rep. (2013). https://CRAN.R-project.org/package=survivalROC.

552

553

40. Tibshirani, R. Regression Shrinkage and Selection via the Lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, 267-288. ISSN: 0035-9246 (1996).

554

555

556

557

41. Eisenhauer, E. A., Therasse, P., et al. New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1). *European Journal of Cancer. Response assessment in solid tumours (RECIST): Version 1.1 and supporting papers* 228-247. ISSN: 0959-8049. doi:10.1016/j.ejca.2008.10.026 (2009).

558

559

560

42. Benzekry, S., Grangeon, M., et al. Machine Learning for Prediction of Immunotherapy Efficacy in Non-Small Cell Lung Cancer from Simple Clinical and Biological Data. *Cancers*, 6210. doi:10.3390/cancers13246210 (2021).

561

562

563

43. Havel, J. J., Chowell, D. & Chan, T. A. The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. *Nature Reviews Cancer*, 133-150. ISSN: 1474-1768. doi:10.1038/s41568-019-0116-x (2019).

564

565

44. Blank, C. U., Haanen, J. B., Ribas, A. & Schumacher, T. N. The "cancer immunogram". *Science*, 658-660. doi:10.1126/science.aaf2834 (2016).

566

567

568

45. Bach, F. R. Bolasso: model consistent Lasso estimation through the bootstrap in (Association for Computing Machinery, New York, NY, USA, 2008), 33-40. ISBN: 978-1-60558-205-4. doi:10.1145/1390156.1390161.

569

570

571

572

46. Peters, S., Dziadziuszko, R., et al. Atezolizumab versus chemotherapy in advanced or metastatic NSCLC with high blood-based tumor mutational burden: primary analysis of BFAST cohort C randomized phase 3 trial. *Nature Medicine*, 1831-1839. ISSN: 1546-170X. doi:10.1038/s41591-022-01933-w (2022).

573

574

575

47. Soyano, A. E., Dholaria, B., et al. Peripheral blood biomarkers correlate with outcomes in advanced non-small cell lung Cancer patients treated with anti-PD-1 antibodies. *Journal for Immunotherapy of Cancer*, 129. ISSN: 2051-1426. doi:10.1186/s40425-018-0447-2 (2018).

576

577

578

579

48. Diem, S., Schmid, S., et al. Neutrophil-to-Lymphocyte ratio (NLR) and Platelet-to-Lymphocyte ratio (PLR) as prognostic markers in patients with non-small cell lung cancer (NSCLC) treated with nivolumab. *Lung Cancer (Amsterdam, Netherlands)*, 176-181. ISSN: 1872-8332. doi:10.1016/j.lungcan.2017.07.024 (2017).

580

581

582

583

49. Peng, L., Wang, Y., et al. Peripheral blood markers predictive of outcome and immune-related adverse events in advanced non-small cell lung cancer treated with PD-1 inhibitors. *Cancer immunology, immunotherapy: CII*, 1813-1822. ISSN: 1432-0851. doi:10.1007/s00262-020-02585-w (2020).

584

585

586

50. Becker, T., Weberpals, J., et al. An enhanced prognostic score for overall survival of patients with cancer derived from a large real-world cohort. *Annals of Oncology*, 1561-1568. ISSN: 0923-7534. doi:10.1016/j.annonc.2020.07.013 (2020).

587

588

589

51. Yang, J.-R., Xu, J.-Y., et al. Post-diagnostic C-reactive protein and albumin predict survival in Chinese patients with non-small cell lung cancer: a prospective cohort study. *Scientific Reports*, 8143. ISSN: 2045-2322. doi:10.1038/s41598-019-44653-x (2019).
590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

52. Bruni, D., Angell, H. K. & Galon, J. The immune contexture and Immunoscore in cancer prog-nosis and therapeutic efficacy. *Nature Reviews Cancer*, 662-680. ISSN: 1474-1768. doi:10.1038/s41568-020-0285-7 (2020).

53. Liu, G., Lu, J., Lim, H. S., Jin, J. Y. & Lu, D. Applying interpretable machine learning workflow to evaluate exposure-response relationships for large-molecule oncology drugs. *CPT: pharmacometrics & systems pharmacology*, 1614-1627. ISSN: 2163-8306. doi:10.1002/psp4.12871 (2022).

54. Desmée, S., Mentré, F., Veyrat-Follet, C. & Guedj, J. Nonlinear Mixed-Effect Models for Prostate-Specific Antigen Kinetics and Link with Survival in the Context of Metastatic Prostate Cancer: a Comparison by Simulation of Two-Stage and Joint Approaches. *The AAPS Journal*, 691-699. ISSN: 1550-7416. doi:10.1208/s12248-015-9745-5 (2015).

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

56. Greillier, L., Monville, F., et al. Abstract LB120: Comprehensive biomarkers analysis to explain resistances to PD1-L1 ICIs: The precision immuno-oncology for advanced non-small cell lung cancer (PIONeeR) trial. *Cancer Research*, LB120. ISSN: 0008-5472. doi:10.1158/1538-7445. AM2022-LB120 (2022).

57. Barlesi, F., Monville, F., et al. Comprehensive biomarkers (BMs) analysis to predict efficacy of PD1-L1 immune checkpoint inhibitors (ICIs) in combination with chemotherapy: a subgroup analysis of the Precision Immuno-Oncology for advanced Non-Small Cell Lung CancER (PIO-NeeR) trial. *Annals of Oncology*. doi:10.1016/iotech/iotech100100 (2022).

58. Assaf, Z. J. F., Zou, W., et al. A longitudinal circulating tumor DNA-based model associated with survival in metastatic non-small-cell lung cancer. *Nature Medicine*, 859-868. ISSN: 1546-170X. doi:10.1038/s41591-023-02226-6 (2023).

59. Nabet, B. Y., Esfahani, M. S., et al. Noninvasive Early Identification of Therapeutic Benefit from Immune Checkpoint Inhibition. *Cell*, 363-376.e13. ISSN: 1097-4172. doi:10.1016/j.cell.2020.09.001 (2020).

60. Cabel, L., Proudhon, C., et al. Clinical potential of circulating tumour DNA in patients receiving anticancer immunotherapy. *Nature Reviews. Clinical Oncology*, 639-650. ISSN: 1759-4782. doi:10.1038/s41571-018-0074-3 (2018).

61. Barrera, L., Montes-Servín, E., et al. Cytokine profile determined by data-mining analysis set into clusters of non-small-cell lung cancer patients according to prognosis. *Annals of Oncology: Official Journal of the European Society for Medical Oncology*, 428-435. ISSN: 1569-8041. doi:10.1093/annonc/mdu549 (2015).

62. Ciccolini, J., Benzekry, S. & Barlesi, F. Deciphering the response and resistance to immune-checkpoint inhibitors in lung cancer with artificial intelligence-based analysis: when PIONeeR meets QUANTIC. *British Journal of Cancer*, 1-2. ISSN: 1532-1827. doi:10.1038/s41416-020-0918-3 (2020).

63. Ciccolini, J., Barbolosi, D., André, N., Barlesi, F. & Benzekry, S. Mechanistic Learning for Combinatorial Strategies With Immuno-oncology Drugs: Can Model-Informed Designs Help Investigators? *JCO Precision Oncology*, 486-491. doi:10.1200/PO.19.00381 (2020).
<table><thead><tr><th>Study</th><th>Description</th><th>Population</th><th>N</th></tr></thead><tbody><tr><td>FIR<br>GO28625</td><td>Phase 2 study for the efficacy and safety of ATZ in advanced NSCLC</td><td>PD-L1 positive locally advanced or metastatic NSCLC (lines 1 and 2+)</td><td>133</td></tr><tr><td>POPLAR<br>GO28753</td><td>Phase 2 randomised controlled trial of ATZ versus docetaxel in NSCLC</td><td>Locally advanced or metastatic NSCLC who failed platinum therapy (line 2)</td><td>134</td></tr><tr><td>BIRCH<br>GO28754</td><td>Phase 2 study of ATZ in advanced or metastatic NSCLC</td><td>Locally advanced or metastatic NSCLC (lines 1, 2 or 3)</td><td>595</td></tr><tr><td>Train</td><td></td><td></td><td>862</td></tr><tr><td>Test - OAK<br>GO28915</td><td>Phase 3 RCT of ATZ versus docetaxel (DTX) in patients with previously treated NSCLC</td><td>Stage IIIB or IV, previously chemo treated</td><td>553</td></tr><tr><td>Train + Test</td><td></td><td></td><td>1415</td></tr></tbody></table>

Four monotherapy studies of atezolizumab in advanced NSCLC. NSCLC: Non-Small Cell Lung Cancer; *p* = number of parameters, *N* = number of patients treated with atezolizumab (patients from French centers were excluded for legal reasons (N=118)). In total, data from 1074 patients from OAK were used as Test set (553 from the ATZ arm, 521 from the DTX arm); PD: Pharmacodynamic; SLD: Sum of the Largest Diameters. CRP: C Reactive Protein; LDH: Lactate Dehydrogenase.

Supplementary Figure 1. Train and test data sets

<table><thead><tr><th>Characteristic</th><th>Total, N = 1415<sup>1</sup></th><th>FIR, N = 133<sup>1</sup></th><th>POPLAR, N = 134<sup>1</sup></th><th>BIRCH, N = 595<sup>1</sup></th><th>OAK, N = 553<sup>1</sup></th><th>p-value<sup>1</sup></th></tr></thead><tbody><tr><td>Age</td><td>64 (57, 70)</td><td>67 (60, 73)</td><td>62 (55, 69)</td><td>65 (57, 71)</td><td>64 (57, 70)</td><td>&lt;0.001</td></tr><tr><td>Sex</td><td></td><td></td><td></td><td></td><td></td><td>0.3</td></tr><tr><td>Female</td><td>568 (40%)</td><td>57 (43%)</td><td>47 (35%)</td><td>251 (42%)</td><td>213 (39%)</td><td></td></tr><tr><td>Male</td><td>847 (60%)</td><td>76 (57%)</td><td>87 (65%)</td><td>344 (58%)</td><td>340 (61%)</td><td></td></tr><tr><td>Weight</td><td>72 (61, 82)</td><td>70 (60, 83)</td><td>73 (63, 84)</td><td>72 (61, 82)</td><td>71 (60, 82)</td><td>0.3</td></tr><tr><td>BMI</td><td>24.9 (22.1, 28.1)</td><td>24.8 (21.9, 27.6)</td><td>25.2 (22.8, 28.7)</td><td>25.0 (22.1, 28.2)</td><td>24.7 (22.0, 28.1)</td><td>0.4</td></tr><tr><td>Unknown</td><td>65</td><td>8</td><td>5</td><td>30</td><td>22</td><td></td></tr><tr><td>Race</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>Asian</td><td>228 (16%)</td><td>6 (4.5%)</td><td>23 (17%)</td><td>77 (13%)</td><td>122 (22%)</td><td></td></tr><tr><td>Others, unknown or missing</td><td>73 (5.2%)</td><td>9 (6.8%)</td><td>9 (6.7%)</td><td>23 (3.9%)</td><td>32 (5.8%)</td><td></td></tr><tr><td>White</td><td>1,114 (79%)</td><td>118 (89%)</td><td>102 (76%)</td><td>495 (83%)</td><td>399 (72%)</td><td></td></tr><tr><td>Smoking history</td><td></td><td></td><td></td><td></td><td></td><td>0.13</td></tr><tr><td>Current</td><td>171 (12%)</td><td>18 (14%)</td><td>22 (16%)</td><td>60 (10%)</td><td>71 (13%)</td><td></td></tr><tr><td>Never</td><td>248 (18%)</td><td>16 (12%)</td><td>27 (20%)</td><td>102 (17%)</td><td>103 (19%)</td><td></td></tr><tr><td>Previous</td><td>996 (70%)</td><td>99 (74%)</td><td>85 (63%)</td><td>433 (73%)</td><td>379 (69%)</td><td></td></tr><tr><td>Characteristic</td><td>Total, N = 1415<sup>1</sup></td><td>FIR, N = 133<sup>1</sup></td><td>POPLAR, N = 134<sup>1</sup></td><td>BIRCH, N = 595<sup>1</sup></td><td>OAK, N = 553<sup>1</sup></td><td>p-value<sup>1</sup></td></tr><tr><td>Heart rate</td><td>81 (71, 92)</td><td>80 (70, 94)</td><td>84 (74, 96)</td><td>80 (70, 90)</td><td>81 (72, 93)</td><td>0.049</td></tr><tr><td>Systolic blood pressure</td><td>122 (111, 133)</td><td>122 (113, 132)</td><td>122 (112, 131)</td><td>120 (110, 133)</td><td>123 (113, 135)</td><td>0.13</td></tr></tbody></table>

Supplementary Figure 2. Patient characteristics: demographics and clinics

<footnote n="1">¹Median (IQR); n (%)</footnote>

<footnote n="2">²Kruskal-Wallis rank sum test; Pearson's Chi-squared test</footnote>
<table><thead><tr><th>Characteristic</th><th>Total, N = 1415<sup>1</sup></th><th>FIR, N = 133<sup>1</sup></th><th>POPLAR, N = 134<sup>1</sup></th><th>BIRCH, N = 595<sup>1</sup></th><th>OAK, N = 553<sup>1</sup></th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>Disease type</td><td></td><td></td><td></td><td></td><td></td><td>0.4</td></tr><tr><td>Locally advanced</td><td>77 (5.4%)</td><td>3 (2.3%)</td><td>8 (6.0%)</td><td>32 (5.4%)</td><td>34 (6.1%)</td><td></td></tr><tr><td>Metastatic</td><td>1,338 (95%)</td><td>130 (98%)</td><td>126 (94%)</td><td>563 (95%)</td><td>519 (94%)</td><td></td></tr><tr><td>Line</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>≥2</td><td>1,255 (89%)</td><td>102 (77%)</td><td>134 (100%)</td><td>466 (78%)</td><td>553 (100%)</td><td></td></tr><tr><td>1</td><td>160 (11%)</td><td>31 (23%)</td><td>0 (0%)</td><td>129 (22%)</td><td>0 (0%)</td><td></td></tr><tr><td>Histology</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>Non-squamous</td><td>1,016 (74%)</td><td>95 (100%)</td><td>87 (65%)</td><td>427 (72%)</td><td>407 (74%)</td><td></td></tr><tr><td>Squamous</td><td>361 (26%)</td><td>0 (0%)</td><td>47 (35%)</td><td>168 (28%)</td><td>146 (26%)</td><td></td></tr><tr><td>Unknown</td><td>38</td><td>38</td><td>0</td><td>0</td><td>0</td><td></td></tr><tr><td>Stage</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>I</td><td>123 (8.9%)</td><td>11 (8.7%)</td><td>4 (3.0%)</td><td>73 (12%)</td><td>35 (6.5%)</td><td></td></tr><tr><td>II</td><td>140 (10%)</td><td>11 (8.7%)</td><td>9 (6.8%)</td><td>73 (12%)</td><td>47 (8.7%)</td><td></td></tr><tr><td>III</td><td>367 (27%)</td><td>30 (24%)</td><td>37 (28%)</td><td>165 (28%)</td><td>135 (25%)</td><td></td></tr><tr><td>IV</td><td>753 (54%)</td><td>75 (59%)</td><td>82 (62%)</td><td>273 (47%)</td><td>323 (60%)</td><td></td></tr><tr><td>Unknown</td><td>32</td><td>6</td><td>2</td><td>11</td><td>13</td><td></td></tr><tr><td>Number of metastases</td><td></td><td></td><td></td><td></td><td></td><td>0.4</td></tr><tr><td>1</td><td>386 (28%)</td><td>35 (27%)</td><td>34 (26%)</td><td>158 (27%)</td><td>159 (29%)</td><td></td></tr><tr><td>2</td><td>655 (48%)</td><td>55 (42%)</td><td>61 (47%)</td><td>289 (50%)</td><td>250 (46%)</td><td></td></tr><tr><td>3</td><td>334 (24%)</td><td>41 (31%)</td><td>34 (26%)</td><td>128 (22%)</td><td>131 (24%)</td><td></td></tr><tr><td>Unknown</td><td>40</td><td>2</td><td>5</td><td>20</td><td>13</td><td></td></tr><tr><td>Liver metastases</td><td>272 (19%)</td><td>29 (22%)</td><td>33 (25%)</td><td>105 (18%)</td><td>105 (19%)</td><td>0.3</td></tr><tr><td>Number of met. loc.</td><td></td><td></td><td></td><td></td><td></td><td>0.3</td></tr><tr><td>Four sites</td><td>70 (4.9%)</td><td>13 (9.8%)</td><td>6 (4.5%)</td><td>27 (4.5%)</td><td>24 (4.3%)</td><td></td></tr><tr><td>One site</td><td>426 (30%)</td><td>37 (28%)</td><td>39 (29%)</td><td>178 (30%)</td><td>172 (31%)</td><td></td></tr><tr><td>Three sites</td><td>264 (19%)</td><td>28 (21%)</td><td>28 (21%)</td><td>101 (17%)</td><td>107 (19%)</td><td></td></tr><tr><td>Two sites</td><td>655 (46%)</td><td>55 (41%)</td><td>61 (46%)</td><td>289 (49%)</td><td>250 (45%)</td><td></td></tr><tr><td>Tumor size</td><td>59 (43, 95)</td><td>59 (57, 59)</td><td>70 (43, 107)</td><td>60 (37, 96)</td><td>70 (43, 102)</td><td>0.001</td></tr></tbody></table>

<table><caption>**Supplementary Figure 3.** Patient characteristics: disease</caption><thead><tr><th>Characteristic</th><th>Total, N = 1415<sup>1</sup></th><th>FIR, N = 133<sup>1</sup></th><th>POPLAR, N = 134<sup>1</sup></th><th>BIRCH, N = 595<sup>1</sup></th><th>OAK, N = 553<sup>1</sup></th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>PD-L1 tumor cells</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>0</td><td>606 (43%)</td><td>4 (3.0%)</td><td>40 (30%)</td><td>212 (36%)</td><td>350 (64%)</td><td></td></tr><tr><td>1</td><td>209 (15%)</td><td>21 (16%)</td><td>56 (42%)</td><td>69 (12%)</td><td>63 (11%)</td><td></td></tr><tr><td>2</td><td>381 (27%)</td><td>105 (79%)</td><td>37 (28%)</td><td>160 (27%)</td><td>79 (14%)</td><td></td></tr><tr><td>3</td><td>217 (15%)</td><td>3 (2.3%)</td><td>1 (0.7%)</td><td>154 (26%)</td><td>59 (11%)</td><td></td></tr><tr><td>Unknown</td><td>2</td><td>0</td><td>0</td><td>0</td><td>2</td><td></td></tr><tr><td>PD-L1 immune cells</td><td></td><td></td><td></td><td></td><td></td><td>&lt;0.001</td></tr><tr><td>0</td><td>252 (18%)</td><td>3 (2.3%)</td><td>40 (30%)</td><td>16 (2.7%)</td><td>193 (35%)</td><td></td></tr><tr><td>1</td><td>398 (28%)</td><td>5 (3.8%)</td><td>56 (42%)</td><td>122 (21%)</td><td>215 (39%)</td><td></td></tr><tr><td>2</td><td>439 (31%)</td><td>25 (19%)</td><td>19 (14%)</td><td>297 (50%)</td><td>98 (18%)</td><td></td></tr><tr><td>3</td><td>322 (23%)</td><td>99 (75%)</td><td>19 (14%)</td><td>159 (27%)</td><td>45 (8.2%)</td><td></td></tr><tr><td>Unknown</td><td>4</td><td>1</td><td>0</td><td>1</td><td>2</td><td></td></tr><tr><td>ECOG</td><td></td><td></td><td></td><td></td><td></td><td>0.6</td></tr><tr><td>Status 0</td><td>505 (36%)</td><td>42 (32%)</td><td>44 (33%)</td><td>215 (36%)</td><td>204 (37%)</td><td></td></tr><tr><td>Status 1 or 2</td><td>909 (64%)</td><td>90 (68%)</td><td>90 (67%)</td><td>380 (64%)</td><td>349 (63%)</td><td></td></tr><tr><td>Unknown</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td></td></tr></tbody></table>

<caption>**Supplementary Figure 4.** Patient characteristics: PD-L1 and ECOG</caption>

<footnote n="1">¹n (%); Median (IQR)</footnote>

<footnote n="2">²Pearson's Chi-squared test; Kruskal-Wallis rank sum test</footnote>

<footnote n="1">¹n (%)</footnote>

<footnote n="2">²Pearson's Chi-squared test</footnote>
Comparative cross-validated predictive performances of the survival ML models.

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

Note: 10-fold cross-validation on FIR, BIRCH and POPLAR (train data set) – performances using all features

15

<caption>**Supplementary Figure 5.** Comparison of ML algorithms and tuning methods</caption>

645

648

<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 6.** TK modeling goodness-of-fit</figcaption></figure>

648

650
<figure><img src="image_9.png" /><figcaption>**Supplementary Figure 7.** Examples of longitudinal kinetics: Albumin</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 8.** Examples of longitudinal kinetics: CRP</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 9.** Examples of longitudinal kinetics: LDH</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 10.** Examples of longitudinal kinetics: Neutrophils</figcaption></figure>
<table><thead><tr><th rowspan="2">Model</th><th colspan="2">Albumin</th><th colspan="2">CRP</th><th colspan="2">LDH</th><th colspan="2">Neutrophils</th></tr><tr><th>BICc</th><th>b</th><th>BICc</th><th>b</th><th>BICc</th><th>b</th><th>BICc</th><th>b</th></tr></thead><tbody><tr><td>Double-exponential</td><td>48,395</td><td>0.058</td><td>28,764</td><td>0.21</td><td>39,886</td><td>0.56</td><td>102,449</td><td>0.14</td></tr><tr><td>Hyperbolic</td><td>48,007</td><td>0.056</td><td>29,712</td><td>0.22</td><td>40,915</td><td>0.62</td><td>102,943</td><td>0.14</td></tr><tr><td>Linear</td><td>49,436</td><td>0.063</td><td>30,020</td><td>0.23</td><td>42,462</td><td>0.70</td><td>105,193</td><td>0.17</td></tr><tr><td>Constant</td><td>49,724</td><td>0.065</td><td>31,332</td><td>0.25</td><td>42,982</td><td>0.74</td><td>108,249</td><td>0.18</td></tr></tbody></table>

Corrected Bayesian Information Criterion (BICc) for four empirical kinetic models of BK. b: standard deviation of the proportional error model

663

665

10

<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 11. Goodness-of-fit metrics of dynamic BK models**</figcaption></figure>

666

668

<caption>**Supplementary Figure 12. Albumin: hyperbolic individual fits**</caption>
Individual residuals

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

⇒ Globally, several individuals not adequately fitted

⇒ Nevertheless, some patients have very accurate fits

⇒ Some trend of the model to underestimation

⇒ Distribution of IWRES remains close to gaussian

Individual fits in patients with at least 10 data points

<figure><img src="image_15.png" /><figcaption>Supplementary Figure 13. CRP: dexp individual fits</figcaption></figure>

Individual residuals

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

⇒ Very good description of individual kinetics

⇒ No sign of model misspecification at
individual level

Individual fits in patients with at least 10 data points

<figure><img src="image_17.png" /><figcaption>Supplementary Figure 14. LDH: dexp individual fits</figcaption></figure>
Individual residuals

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

⇒ Moderate description of individual kinetics

⇒ No major sign of model misspecification at
individual level

Individual fits in patients with at least 10 data points

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

<figure><img src="image_20.png" /><figcaption>**Supplementary Figure 15.** Neutrophils: dexp individual fits</figcaption></figure>

<caption>**Supplementary Figure 16.** ROC curves for variables landmark times (test set - OAK)</caption>
<figure><img src="image_21.png" /><figcaption>**Supplementary Figure 17.** Learning curve</figcaption></figure>

17

681

683

<figure><img src="image_22.png" /><figcaption>**Supplementary Figure 18.** Additional value of NLME to baseline for multiple metrics</figcaption></figure>

684

685
<figure><img src="image_23.png" /></figure>

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

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

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

Feature set
* BK
* mBBL
* TX

<caption>**Supplementary Figure 19.** kML models using only single kinetic markers</caption>
691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

## Dimensionality reduction for RNAseq

Initial expression data from RNAseq consisted of 715 patients and 58,311 transcripts. The first step of data filtering removed all transcripts with less than 10 read counts for all patients, then selected genes with highest variability between patients (top 15,000 transcripts most variable). Then, data were normalized using upper quartile normalization which consisted in dividing each read count by the 75th percentile of the read counts of the corresponding sample and the final expression values were log₂ transformed. Subsequently, a univariable Cox regression model was employed to statistically assess the correlations between the expression levels of the transcripts and overall survival. Bonferroni correction was used to adjust p-values from multiple univariate tests. This step was performed using the RegParallel R package. We selected transcripts with high predictive values using following criteria: adjusted log rank < 0.01 and HR < 0.85 or HR > 1.2. The remaining transcripts were used to perform a bootstrap Lasso Cox regression with cross-validation using mainly the glmmnet R package. Finally, the smallest number of transcripts with best predictive model (highest C-index) was selected for further analysis.

## Rules for BK processing

1. Observations outside lower (LB) and upper (UB) physiological bounds were discarded using the following values, determined from discussion with a clinical oncologist: albumin, LB = 10 g L⁻¹, UB = 100 g L⁻¹; CRP, no LB, UB = 300 mg L⁻¹; LDH, LB = 50 U/L, UB = 2000 U/L; neutrophils, no LB, UB = 20.

2. For duplicates, the first one recorded was kept.

3. Denoting BKₙ the value of the a "BK" at time tₙ for a given patient, we excluded values such that: BKₙ ∉ (BKₙ⁻₁, BKₙ₊₁) AND |BKₙ − BKₙ⁻₁| > 3 × sd<sub>BK</sub> AND |BKₙ − BKₙ₊₁| > 3 × sd<sub>BK</sub>, where sd<sub>BK</sub> is the standard deviation of {BKₙ}ₙ, i.e. all time points for this patient.

4. The BK value at the closest time point to treatment initiation was kept, provided this time point was no more than 40 days before or 10 days after treatment initiation (otherwise, patient was disregarded).

## Nonlinear mixed-effects modeling

Denoting by M(t; θ) a structural dynamic model that depends on time t and a set of parameters θ, longitudinal observations y<sup>i</sup><sub>j</sub> in patient i at time t<sup>i</sup><sub>j</sub> were assumed to follow the observation model

$$y_j^i = M(t_j^i; \theta^i) + \epsilon_j^i,$$

where ε<sup>i</sup><sub>j</sub> ~ N(0, σ<sup>i</sup><sub>j</sub>) is the gaussian-distributed error model. The latter was either constant (σ<sup>i</sup><sub>j</sub> = a, ∀i, j) for TK or proportional (σ<sup>i</sup><sub>j</sub> = bM(t<sup>i</sup><sub>j</sub>; θ<sup>i</sup>), ∀i, j) for BK. To describe inter-individual variability, individual parameters θ<sup>i</sup> were assumed to follow log-normal distributions:

$$\ln(\theta^i) = \ln(\theta_{pop}^i) + \eta^i, \quad \eta^i \sim N(0, \omega^2)$$

with population-level parameters θ<sub>pop</sub> and ω. Estimation of these was performed using the stochastic approximation of expectation maximization algorithm implemented in the Monolix software.