M Kerioui
email: marion.kerioui@inserm.fr

S Desmée

J Bertrand

F Mercier

J Y Jin

R Bruno

J Guedj

Nonlinear multilevel joint model for individual lesion kinetics and survival to characterize intra-individual heterogeneity in patients with advanced cancer

Keywords: Joint model, Multilevel modelling, Nonlinear mixed-effects model, Bayesian inference, Advanced cancer

In advanced cancer patients, tumor burden assessment relies on the Sum of the Longest Diameters (SLD) of the target lesions, a marker that lumps all lesions together and ignores intra-patient heterogeneity. Here, we relied on a rich dataset of 342 metastatic bladder cancer patients treated with a novel immunotherapy agent to develop a Bayesian multilevel joint model that can quantify the heterogeneity in lesion dynamics and measure their impact on survival.

Using a nonlinear model of tumor growth inhibition, we estimated that dynamics differed greatly among lesions, and inter-lesion variability accounted for about 35% of the total variance of both tumor shrinkage and treatment effect duration. Next, we investigated the impact of individual lesion dynamics on survival. Lesions located in the liver and in the bladder had twice as much impact on the instantaneous risk of death as compared to those located in the lung or the lymph nodes. Finally we evaluated the gain of individual lesion follow-up for dynamic predictions.

Consistent with results at the population levels, the individual lesion model outperformed a model relying only on SLD, especially at early landmark times and in patients having liver or bladder target lesions. Our results show that the use of SLD leads to a loss of information and our model can be used to characterize tumor dynamics and survival of advanced cancer patients.

Introduction

The measure of tumor burden relies on the Sum of the Longest Diameters (SLD) of the target lesions, following the RECIST criteria [START_REF] Eisenhauer | New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1)[END_REF]. In advanced cancer patients, where several target lesions, located in distinct organs, may be present, the use of SLD therefore severely reduces the information available by lumping together all measured lesions. This drawback, which has long been identified in the oncology and statistical literature [START_REF] Dong | Mixed Responses to Systemic Therapy Revealed Potential Genetic Heterogeneity and Poor Survival in Patients with Non-Small Cell Lung Cancer[END_REF], has become more problematic with the advent of immunotherapy treatments, which could be associated with an increased heterogeneity in the individual lesions to treatment (Tozuka et al., 2020;[START_REF] Humbert | Dissociated Response in Metastatic Cancer: An Atypical Pattern Brought Into the Spotlight With Immunotherapy[END_REF].

In the past decade, semi-mechanistic nonlinear joint models of tumor kinetics and survival have been proposed to characterize the response to treatment and its impact on survival [START_REF] Desmée | Using the SAEM algorithm for mechanistic joint models characterizing the relationship between nonlinear PSA kinetics and survival in prostate cancer patients: Joint Model for Nonlinear Kinetics and Survival Data[END_REF][START_REF] Król | Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death: Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death[END_REF][START_REF] Gavrilov | 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[END_REF][START_REF] Tardivon | Association Between Tumor Size Kinetics and Survival in Patients With Urothelial Carcinoma Treated With Atezolizumab: Implication for Patient Follow-Up[END_REF], while limiting the bias caused by informative censoring [START_REF] Desmée | 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[END_REF][START_REF] Mbogning | Joint modelling of longitudinal and repeated time-to-event data using nonlinear mixed-effects models and the stochastic approximation expectation-maximization algorithm[END_REF][START_REF] Björnsson | Performance of nonlinear mixed effects models in the presence of informative dropout[END_REF]. Performant algorithms have been developed, relying on softwares such as Monolix or NONMEM, allowing the diffusion of these methods in the statistical community (Kerioui et al., 2022). Joint models accounting for individual lesion dynamics have recently been developed, using existing softwares such as the R package rstanarm, but this approach remains limited to linear models for the longitudinal parts [START_REF] Brilleman | Joint longitudinal and timeto-event models for multilevel hierarchical data[END_REF]. Conversely, nonlinear models have been developed for individual lesion dynamics, but these models did not integrate the association with survival [START_REF] Mercier | Longitudinal analysis of organ-specific tumor lesion sizes in metastatic colorectal cancer patients receiving first line standard chemotherapy in combination with anti-angiogenic treatment[END_REF] or remained limited to a two-stage inference approach [START_REF] Krishnan | Tumor growth inhibition modeling of individual lesion dynamics and interorgan variability in HER2-negative breast cancer patients treated with docetaxel[END_REF]. Developing nonlinear joint models for individual lesions is made complex by the fact that no existing softwares can readily be extended to incorporate another layer of random effects capturing intra-patient variability. In a frequentist approach, the curse of dimensionality may be a main limitation to its implementation. Beside the numerical complexity, using individual January 2022 lesions requires to have access to large dataset, where individual lesions are scrutinized and reported in a large number of patients. Here, we built on a previous works showing that the Hamiltonian Monte-Carlo No U-Turn Sampler (HMC-NUTS algorithm) implemented in Stan software could provide reliable estimates of nonlinear joint models [START_REF] Kerioui | Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy[END_REF]. We extended this approach by taking into account the variability across lesions within a same individual, considering that this variability was caused by both different organ location as well as an intrinsic within-patient variability. We also relied on a rich data of data from 342 advanced bladder cancer patients treated with immune-checkpoint inhibitor agent atezolizumab [START_REF] Powles | Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial[END_REF], and for which up to 5 lesions were followed over time (Kerioui et al., 2022). We showed that this model could well characterize the variability in individual lesion dynamics, and we evaluated its role to improve survival's prediction, relying on individual dynamic approaches (Desmée et al., 2017a).

Methods

Data

IMvigor211 is a phase 3 randomized clinical trial, including second-line patients suffering from advanced or metastatic bladder cancer and treated with anti PD-L1 monoclonal antidbody atezolizumab, versus a chemotherapy control arm [START_REF] Powles | Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial[END_REF].

A total of 459 patients were treated with atezolizumab, and 2 patients were excluded of our analysis as they did not have any target lesion assessment. We focused on patients having only target lesions in the four main organs, namely the lymph nodes (1), the lung (2), the liver (3) and the bladder (4), leading to a final analysis population of N = 342 patients.

For each patient, longitudinal follow-up of individual target lesions and survival data were available.

Tumor structural model

We used the simplified Tumor Growth Inhibition (sTGI) model [START_REF] Claret | Evaluation of Tumor-Size Response Metrics to Predict Overall Survival in Western and Chinese Patients With First-Line Metastatic Colorectal Cancer[END_REF] to describe T S(t) the tumor size at time t. At patient inclusion t = 0, the tumor size is T S 0 and the model assumes an exponential growth g of the tumor size in absence of treatment.

At treatment initiation (t = t x ), the model assumes a tumor growth inhibition ǫ, and the treatment effect decreases over time (i.e. due to treatment resistance appearance) with rate c, as followed:

T S(t) =          T S 0 × e g×t t < t x T S 0 × exp g × t -ǫ c × (1 -e -c×(t-tx) ) t t x
(1)

Nonlinear mixed-effects model

Let y i,j,k,l denote the l th measurement of the k th target lesion in organ j ∈ {1, . . . , 4} in patient i ∈ {1, . . . , N }. We assume a limit of quantification loq = 2.5mm below which the lesion cannot be properly measured. We denote K i,j the number of target lesions followed over time in organ j of patient i, and K i,• = 4 j=1 K i,j the total number of target lesions for patient i. We denote L i,j,k the number of measurements for lesion k in organ j of patient i.

The mixed-effect model was defined as:

y i,j,k,l = T S(t i,j,k,l , ψ i,j,k ) + (σ 1,j + σ 2,j × T S(t i,j,k,l , ψ i,j,k )) e i,j,k,l (2)
with T S(t i,j,k,l , ψ i,j,k ) the structural model function defined in equation 1 that describes measurement y i,j,k,l at time t i,j,k,l , depending on ψ i,j,k = {ψ i,j,k,T S 0 , ψ i,j,k,ǫ , ψ i,j,k,g , ψ i,j,k,c } the individual parameters vector specific to the lesion k in organ j of patient i. The parameters σ 1,j and σ 2,j are respectively the additive and multiplicative components of the error term specific to organ j, and e i,j,k,l is a residual error term assumed to follow a centered gaussian January 2022 distribution.

The individual parameters are assumed to follow a log-normal distribution, as followed:

log(ψ i,j,k ) = log(µ) + ξ j + η i + ρ i,j,k (3)
where µ is a fixed-effect parameters vector and ξ j is an organ fixed-effects vector specific to organ j. To ensure identifiability of the parameters, a null sum constraint was assumed on the organ fixed-effect, i.e. 4 j=1 ξ j = 0.

The first random effect is specific to patient i:

η i ∼ N (0, Ω 1 )
and Ω 1 is a diagonal variancecovariance matrix characterizing the inter-patient variability. The second random effect is the lesion effect: ρ i,j,k ∼ N (0, Ω 2 ) and Ω 2 is a diagonal variance-covariance matrix characterizing the inter-lesions variability. For all patients i and all lesions k of all organs j, we assume independence between η i and ρ i,j,k .

The proportion of total variance accounted for by the patient and the lesion level can be derived for each biological parameter as the diagonal components of the matrix

Ω 1 (Ω 1 + Ω 2 ) -1
and Ω 2 (Ω 1 + Ω 2 ) -1 respectively.

Joint model

We denote (T i , δ i ) the vector of survival data of patient i, with T i the time of event for patient i, corresponding to the minimum between the time of death X i and the censoring time C i , and

δ i = 1 {X i C i } the indicator of death. The individual hazard function h i (t) describes the survival data of patient i (T i , δ i ): h i (t|θ, ψ i ) = h 0 (t) exp (β × f (t, ψ i )) (4)
with h 0 (t) = γ λ ( t λ ) γ-1 a Weibull baseline hazard function with shape parameter γ and scale parameter λ, β the vector of link parameters, which characterizes the strength of the association between the longitudinal biomarkers and the survival process, and f (t, ψ i ) the vector of link functions, which depends on time and the vector of individual parameters specific to patient i, ψ i = (ψ i,j,k ) j∈{1,...,4}, k∈{1,...,K i,j } . In the following, we denote θ = {µ, ξ, Ω 1 , Ω 2 , σ 1 , σ 2 , γ, λ, β} the vector of population parameters for the joint model. Of note, the individual survival function S i (t|θ, ψ i ) can be directly derived from the individual hazard function as: 1 gives the Directed Acyclic Graph (DAG) of the joint model.

S i (t|θ, ψ i ) = exp - t 0 h i (u|θ, ψ i )du . Figure
[Figure 1 about here.]

Association between tumor size and survival

We considered different association structures between the individual lesions and the instantaneous risk of death. The first link function assumes the same impact β lesion of each target lesion on the hazard function, and is called the "tumor burden" model, defined as follows:

β lesion × 4 j=1 K i,j k=1 T S(t, ψ i,j,k ) (5)
The second link function considers an interaction between organ and tumor size and is called the "organ tumor burden" model, defined as follows:

4 j=1 β j × K i,j k=1 T S(t, ψ i,j,k ) (6)
The third link function considers the impact of the maximum size of the lesions and is called the "maximum" model, defined as follows:

β max × max j,k {T S(t, ψ i,j,k )} (7)
where max j,k T S(t, ψ i,j,k ) is the size of the largest target lesion in patient i and time t.

January 2022

Finally, the fourth link function incorporates the impact of intra-patient heterogeneity, using the difference between the maximum and the minimum of the lesions at each time, and is called the "range" model, defined as follows:

β range × max j,k {T S(t, ψ i,j,k )} -min j,k {T S(t, ψ i,j,k )} (8)
where min j,k T S(t, ψ i,j,k ) is the size of the smallest target lesion in patient i at time t.

Model estimation

Conditionally to the population parameters and the random effects, we assume independence between all longitudinal measurements of one patient, as well as independence between the longitudinal and survival data of one patient. Hence, denoting y i the vector of all lesions measurements in patient i, the individual likelihood is the product of the densities and can be defined as follows:

L(y i , T i , δ i |θ, ψ i ) = 4 j=1 K i,j k=1 L i,j,k l=1 p L (y i,j,k,l |θ, ψ i,j,k ) × p S (T i , δ i |θ, ψ i ) (9)
where p L (y i,j,k,l |θ, ψ i,j,k ) is the density of a normal distribution with mean T S(t i,j,k,l , ψ i,j,k )

and variance σ 1,j + σ 2,j × T S(t i,j,k,l , ψ i,j,k ), evaluated in y i,j,k,l , and p S (T i , δ i |θ, ψ i ) is the density of a censored survival process with hazard function h i (t|θ, ψ i ), defined as follows:

p S (T i , δ i |θ, ψ i ) = h δ i i (T i |θ, ψ i ) × S i (T i |θ, ψ i ) (10)
Of note, the lesion measurements below the limit of quantification were considered as left-censored in the longitudinal density [START_REF] Jacqmin-Gadda | Analysis of leftcensored longitudinal data with application to viral load in HIV infection[END_REF] (see section 10 in supplementary material for more details).

Given the vector of population parameters θ, the random effects of any two lesions and the random effects of any two patients are independent. Denoting η and ρ the vector of the patients and lesions random effects respectively, we can express the prior distribution on the parameters as:

p(θ, η, ρ) ∝ N i=1 4 j=1 K i,j k=1 p(ρ i,j,k |θ)p(η i |θ)p(θ) (11)
Here, the prior distributions of the random effects are given by the model specification, i.e.

η i |θ ∼ N (0, Ω 1 ) and ρ i,j,k |θ ∼ N (0, Ω 2 ).
According to Bayes' theorem, the posterior distribution can be defined as follows:

p (θ, η, ρ|y, T , δ) ∝ N i=1 4 j=1 K i,j k=1 L i,j,k l=1 p L (y i,j,k,l |θ, ψ i,j,k )p S (T i , δ i |θ, ψ i )p(ρ i,j,k |θ)p(η i |θ)p(θ)
All parameters were estimated using HMC algorithm in Stan via the R package rstan (version 2.19.2). We ran three chains in parallel of 500 iterations each, including 250 warmup iterations, leading to a posterior sample of P = 750 replicates. To compute the survival function, the integral of the hazard function was approximated by a Gauss-Kronrod quadrature method. Stan code is available in supplementary material.

Convergence was assessed based on the observation of posterior trace plots, and classic Stan diagnostic tools, including R and effective sample sizes. The classic split-R compares the within and between-chain variances, and is expected to be close to 1 to ensure good convergence of the chains. The effective sample size is a measure of sample effectiveness, the larger the better.

Prior distributions

Based on previous analyses and parameters biological definitions, we set informative prior distributions on the fixed-effect parameters µ, i.e. µ T S 0 ∼ N (25, 5), µ ǫ ∼ Γ(1, 100), µ g ∼ January 2022

Γ(1, 100) and µ c ∼ Γ(1, 1000). All random effects variance priors were half-Cauchy distributions, with Cauchy(0, 1) + for diagonal terms of Ω 1 and Cauchy(0, 0.5) + for diagonal terms of Ω 2 . Other prior distributions were weakly informative, i.e. ξ ∼ N (0, I), with I the identity matrix, σ 1,j ∼ LN (0, 1), σ 2,j ∼ LN (0, 1), γ ∼ LN (0, 1), λ ∼ N (500, 500) + and all components of link parameters vector β follows a N (0, 1000).

Model comparison

For each association structure between the individual lesions and the survival, the Bayesian criteria WAIC was computed. The better the model fits to the data, the lowest this criterion is. We defined the WAIC as follows:

W AIC = -2 × N i=1 log E post (L(y i , T i , δ i |θ, ψ i )) + 2 × p W AIC (12) p W AIC = N i=1 V ar post (log L(y i , T i , δ i |θ, ψ i )) (13)
where

E post (L(y i , T i , δ i |θ, ψ i ))
is the posterior mean of the individual likelihood and can be estimated as:

E post (L(y i , T i , δ i |θ, ψ i )) = 1 P P p=1 L(y i , T i , δ i |θ p , ψ p i ), (14)
and V ar post (log L(y i , T i , δ i |θ, ψ i )) is the posterior variance of the individual log-likelihood and can be estimated as:

V ar post (log L(y i , T i , δ i |θ, ψ i )) = V P p=1 (log L(y i , T i , δ i |θ p , ψ p i )), (15)
where V P p=1 represents the sample variance.

Posterior Predictive Checks

Posterior predictive checking (PPC) relies on the visual comparison of the original dataset with many replicated data generated under the selected model.

We replicated B = 1000 datasets using the following procedure. For each replicated dataset b, we drew a vector of population parameters in the posterior distribution θ b ∼ p(θ|y, T , δ)

and we sampled the random effects of each replicated lesion k from organ j of each replicated patient i in their distributions η b i ∼ N (0, Ω b 1 ) and ρ b i,j,k ∼ N (0, Ω b 2 ). The individual parameters vector of patient i ψ i b was then computed as function of θ b , η b i and ρ b i,j,k following equation ( 3). Following study protocol, for each theoretical measurement time t l , we derived the simulated trajectory of the lesion, drawing

y b i,j,k,l ∼ p L (y i,j,k,l |θ b , ψ b i,j,k ),
and we sampled a time-to-death

T b i ∼ p S (T i |θ b , ψ b i ).
In order to be consistent with study protocol, we assumed that SLD was no longer measured in case of progression defined by an increase of more than 20% and more than 5 mm from nadir. We also kept the same structure than the observed data, in terms of number of lesions per patient and location of the lesions.

For each time t l , and depending on the lesion location, we computed the median of the lesions measurements over the N simulated patients. We used the 2.5 th and 97.5 th percentiles of the medians over the B replicated datasets to provide 95% prediction intervals and we compared this prediction with the median of the lesion measurements in the original dataset. Similarly, we compared the survival probability curve estimated by Kaplan-Meier on the original dataset with the 95% prediction interval computed over the B Kaplan-Meier estimates of the survival probability in each B dataset.

Reference joint model of SLD and survival

We compared the predictive performances of our final model with a model that ignores lesion heterogeneity and simply lumps all lesions together [START_REF] Kerioui | Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy[END_REF]) (comparison criteria are described in the following, section 2.12). In this reference model, only the Sum of the Longest Diameters (SLD) of the target lesions was fitted for the longitudinal part, and we assumed an association between the current value of the SLD and the instantaneous risk of death (see section 4 in supplementary material for a formal definition of the reference model).

January 2022

Individual dynamic predictions

Finally, we aimed to assess the model ability to predict the survival probability of a new patient, given his longitudinal history up to a certain time.

We assumed that the vector of population parameters θ was known. We were interested in the ability of the model to predict the future of a patient i still at risk of death at a landmark time s, based on his longitudinal follow-up until s, Y i (s) = {y i,j,k,l ; 0 t i,j,k,l s}. More specifically, we aimed to predict the conditional probability of death between s and s + t for patient i:

π i (s + t|s) = P (X i < s + t|X i > s, Y i (s), θ) (16)
where t > 0 is called the horizon time. To this purpose, at each landmark time s, we sample M realizations from the conditional distribution of the random effects η i and ρ i,j,k , i.e.:

p(η i , ρ i |X i > s, Y i (s), θ) ∝ 4 j=1 K i,j k=1 L i,j,k (s) l=1 p L (y i,j,k,l |θ, η i , ρ i,j,k )S i (s|θ, ψ i )p(ρ i,j,k |θ)p(η i |θ). ( 17
)
where ρ i = (ρ i,j,k ) j∈{1,...,4},k∈{1,...,K i,j } and L i,j,k (s) is the number of measurements of lesion k in organ j in patient i observed up to time s. In practice, M was set to 600. Estimates of the conditional probabilities of death π i (s + t|s) were then derived (Desmée et al., 2017a) for landmark times s = 0, 90 and 180 days.

Model calibration and discrimination

Discrimination was assessed with the time-dependent Area Under the ROC Curve (AUC) defined as follows:

AU C(s, t) = P (π i (s + t|s) > π j (s + t|s)|s < X i < s + t, X j > s + t) (18)
Calibration was assessed using the time-dependent Brier Score (BS), defined as follows:

BS(s, t) = E (1 s<X<s+t -π(s + t|s)) 2 |X > s . ( 19
)
The larger the AUC, the better the model discriminates between patients at risk of death.

Conversely, the smaller the BS, the better the model is calibrated. The Inverse Probability of Censoring Weighting approach was applied to compute AUC and BS, accounting for censoring bias [START_REF] Blanche | Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks[END_REF][START_REF] Blanche | Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks: Comparing Dynamic Predictive Accuracy of Joint Models[END_REF]. Time-dependent AUC was computed using the timeROC R package, and time-dependent BS with the BS function available in [START_REF] Blanche | Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks: Comparing Dynamic Predictive Accuracy of Joint Models[END_REF].

Results

Data

A total of 751 target lesions were followed over time, including 296 lesions in the lymph, 205 in the lung, 168 in the liver and 82 in the bladder. Each patient had 2 (1-5) target lesions in median (minimum-maximum). A total of 2726 target lesions measurements was collected, with a median (minimum-maximum) of 6 measurements per patient (1-35). Overall, 226 patients (66%) died during the follow-up. Of note, the proportion of death was larger in the 159 patients having at least one target lesions in the liver or the bladder (77%) as compared to the 183 patients without any liver or bladder target lesions (57%).

Intra-patient variability

As results were consistent across the candidate models (see Tables S2,S3 and S4 in supplementary material), and the organ-tumor burden model provided the best fit to the data (see below section 3.4), we mainly focused on this latter in this section and the following.

Individual fits showed good adjustment of the model to the data, despite the large diversity of longitudinal profiles. More specifically, the model was able to capture different lesions kinetics within one patient, including lesions with heterogeneous responses to treatment, i.e. some lesions growing and some other shrinking within one organ of one patient (e.g. patient 1304, in Figure 2).

[Figure 2 about here.]

After accounting for organ effects and inter-patient variability (IPV), the inter-lesion variability (ILV) was non-negligible.

Notably, the ILV on T S 0 (95% CI: [0.33;0.39]) was greater than the IPV (95% CI: [0.33;0.39]), meaning that the heterogeneity was greater between the baseline lesions size of one same patient than between the baseline lesions size of any two patients. Therefore, for the baseline tumor size, the proportion of variance explained by the lesion level was equal to 59%.

The IPV and ILV was particularly large on the treatment related parameters ǫ (IPV 95%

CI: [1.05;1.58] and ILV 95% CI: [0.53;0.81]) and c (IPV 95% CI: [0.79;1.92] and ILV 95% CI:

[0.26;1.24]), which characterize respectively the strength and the durability of the treatment effect. For these parameters, the proportion of total variance explained by the ILV was about 35%.

Finally, the ILV on the g parameter, that characterizes the tumor growth, was more modest (95% CI: [0.09;0.47]), and only accounted for 26% of the total variance, suggesting that the IPV (95% CI: [0.63;1.02]) was to predominate at the later follow-up times (see Table 1).

Organ-specific lesions kinetics

Table 1 gives the posterior mean estimates and 95% CI of the longitudinal population parameters for the organ tumor burden model.

The baseline tumor lesion size was significantly smaller in the lymph (ξ T S 0 ,lymph 95% CI: [-0.30,-0.18]) and in the lung (ξ T S 0 ,lung 95% CI: [-0.24,-0.12]) than in any other location, leading to a baseline lesion size of 21.3 mm and 22.5 mm respectively. Conversely, the baseline tumor lesion size was significantly greater in the bladder (ξ T S 0 ,bladder 95% CI: [0.32,0.48]) than in any other organ, leading to a baseline lesion size twice as large in the bladder as in the lymph nodes (40.3 mm).

The decrease of the lesion size induced by the treatment effect was significantly faster in the lymph nodes (ξ ǫ,lymph 95% CI: [0.16;0.70]) as compared to any other organ., with a lesion shrinkage rate of 0.0038 day -1 . Consequently, we observed an important decrease of the lymph nodes size after treatment initiation (see Figure 3).

Conversely, the lesions shrinkage was more modest in the liver lesions with a shrinkage rate of 0.0018 day -1 , but the duration of the treatment effect was significantly longer in this organ (ξ c,liver 95% CI: [-2.43;-0.68]), for lesions responding to treatment.

Finally, the lesion growth was significantly slower in the lymph (ξ g,lymph 95% CI: [-0.61,-0.21])

and significantly faster in the liver lesions (ξ g,liver 95% CI: [0.38;0.78]) leading to a natural tumor growth twice as fast in the liver (0.0055 day -1 ) as in the lymph (0.0021 day -1 ).

[Table 1 about here.]

[Figure 3 about here.]

Association between individual lesions and survival

Table 2 gives the WAIC criterion of each candidate models and the posterior mean and 95% CI of the link parameters.

[Table 2 about here.]

Under the tumor burden model, the current value of each individual lesion was significantly associated with the instantaneous risk of death, and more specifically a larger lesion size induced a poorer survival probability (β lesion = 0.011 95% CI: [0.0090;0.013]).

In the organ tumor burden model, the impact of the liver lesions on the instantaneous risk of death was larger than in any other organ (β liver 95% CI: [0.011;0.016]), while the impact Both longitudinal and survival PPC showed the ability of the model to correctly replicate the data (Figure 3). Regarding the longitudinal data, the PPC was stratified by organ, and therefore showed the model ability to precisely describe the tumor organ effect on the lesion kinetics. Following previous work from (Kerioui et al., 2022), patients with at least one target lesion in the liver or the bladder were expected to have a poorer survival probability. We therefore stratified the survival PPC by liver/bladder lesion status. Therefore, the model was able to replicate the survival curve in the total population (see Figure S1 in supplementary material), and particularly showing its ability to discriminate survival depending on tumor location (Figure 3).

Individual dynamic predictions

The tumor burden, tumor burden and range and organ tumor burden models gave similar predictive performances in terms of discrimination and calibration (see Figure S3 in supplementary material), although we slightly improved the predictions when including a marker of intra-patient heterogeneity in the link function. As an illustration, at a landmark time s = 90 days and a horizon time t = 90 months, AUC was about 0.73 under the tumor burden model, and 0.76 under the organ tumor burden model. Overall, we observed a better discrimination and calibration (see Brier Scores in Figure S2 in supplementary material) using the organ tumor burden model in comparison to the joint model only relying on the SLD, at the early landmarks s = 0 and 90 days (see Figure 4). At comparison to the AUC under the SLD model at horizon times of 60 and 120 days.

These differences were also slightly exacerbated in the subpopulation of patients having more than one target lesion followed over time (see Figure S4 in Supplementary material).

[Figure 4 about here.]

Discussion

We here developed a multilevel nonlinear joint model for individual lesions kinetics and survival, and inference was conducted in a Bayesian framework using the HMC-NUTS algorithm in Stan software. The model could capture the large variability in individual lesion dynamics, caused by both different organ locations and intrinsic variability. Parameters such as the duration of the treatment effect could vary by more than 10-fold across the organs, and the shrinkage of the tumor size induced by the treatment was as twice as larger in the lymph nodes as compared to the organs with a poorer response, such as the liver. After adjusting for organ location, we estimated that intra-patient variability accounted for about a third of the total parameter variability in the population. Individual lesions also had a differential effect on survival, with lesions located in the liver and in the bladder having twice as much impact on the instantaneous risk of death as compared to those located in the lung or the lymph nodes. At the individual level, our model outperformed a model relying only on SLD and not on the individual lesions, in particular at early landmark times and in patients having target lesions in the liver or the bladder.

Inference was associated with intensive computational duration, and was largely impacted by the complexity of the link function considered. Here, we only used the between-chain parallelization ; however recent advances in Stan to allow the within-chain parallelization, could drastically reduce computation times. Beside calculation time, one major limitation of this model is that it only relies on target lesions, which are defined by the clinician at the patient inclusion, and limited to a maximum of 5 lesions and 2 lesions in one organ per RECIST1.1. This induces a loss of information, which could be addressed by including more lesions in the analysis, but at the expenses of more complex follow-up visits. How many lesions are needed to characterize intra-patient's variability is unknown, but this model could be useful to go beyond the rules of RECIST1.1. In that perspective, this type of approach would also benefit from integrating the effect of new lesions [START_REF] Król | Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death: Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death[END_REF]Zecchin et al., 2016).

To conclude, this model could be used to better understand the sources of variability in the response to treatments. It can be used to follow individual lesion kinetics and better predict their impact on survival, at both population and individual levels.

Table 1

Estimates of population parameters of the organ tumor burden model, posterior mean [95% credibility interval] of the fixed-effects in each organ, inter-patient variability (IPV) and inter-lesion variability (ILV), for each biological parameter, namely the baseline tumor size (T S0), the treatment induced shrinkage (ǫ), the natural tumor growth (g), and the resistance appearance (c). The fixed-effects specific to each organ were computed as functions of the general fixed-effects µ and the organ effects ξj, as defined in equation ( 3).

Fixed-effect

Comparison of candidate model parameter estimates

Tables S2 andS3 show the estimates and 95% CI of longitudinal parameters of the tumor burden model and the tumor burden and range model respectively, compared to the estimates of the organ tumor burden model parameters (Table S4). performances were similar among the three candidate although we observed a slight improvement when including a marker of intra-patient heterogeneity, as compared to the tumor burden model.

Time-dependent AUC and Brier Scores in patients having more than one target lesion

We compared the time-dependent AUC and Brier Scores computed in the total analysis population versus in the subpopulation of patients having more than one target lesion (Figure S4). In this subpopulation, the discreprancies between the individual lesions and SLD follow-up were slightly exacerbated, especially for early horizon times. As an illustration, in the total population and at a landmark time s = 3 months, the organ tumor burden model provided AUC of 0.73 for an horizon time of 5 months and up to 0.76 for an horizon time of 1 month. In comparison, the SLD model provided AUC of 0.70 and 0.73 respectively. In comparison, in the subpopulation, the AUC were about 0.84 and 0.81 under the organ tumor burden model and only 0.74 and 0.77 under the SLD model. We note that the AUC and Brier Scores were better overall when focusing on this subpopulation, suggesting that it is easier to anticipate the future of patients with many target lesions.

Figure S4: Time-dependent AUC (left panels) and Scores (right panels) in the subpopulation of patients having more than one target lesion followed over time, for each landmark time; s = 0 months, s = 3 months, s = 6 months, s = 9 months, and s = 12 months (from lightest to darkest respectively), for the organ tumor burden model (solid lines) or SLD model (dashed lines).

Accounting for the limit of quantification

The lesion measurements below the limit of quantification loq were considered as left-censored in the longitudinal density, and we used the method developed by Jacqmin-Gadda et al. to take into account their contribution to the likelihood. Thus, we denoted ∆ i,j,k,l the left-censor indicator that is 1 if the l th measurement of lesion k in location j in patient i is left-censored, 0 if it is not, and the density of longitudinal data p L is defined as follows:

p L (y i,j,k,l |θ, ψ i,j,k ) = g i,j,k,l (y i,j,k,l |θ, ψ i,j,k ) (1 -∆ i,j ) + Φ i,j,k,l (loq|θ, ψ i,j,k )∆ i,j,k,l

where g i,j,k,l is the probability density and Φ i,j,k,l the distribution function of a normal distribution with expectation T S(t i,j,k,l , ψ i,j,k ) and variance σ 1,j + σ 2,j × T S(t i,j,k,l , ψ i,j,k ), evaluated in y i,j,k,l .

a

  landmark time s = 90 days, the organ tumor burden model provided AUC of 0.73 for an horizon time of 150 days and up to 0.76 for an horizon time of 30 days. In comparison, the SLD model provided AUC of 0.70 and 0.73 respectively. The differences between the organ tumor burden model and the SLD model were exacerbated in the subpopulation of patients having at least one liver or bladder target lesion. As an illustration, at a landmark time s = 90 days, the AUC under the organ tumor burden model was improved by 0.10 points in January 2022

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Directed Acyclic Graph (DAG) of the multi-hierarchical joint model, including population (yellow), patient (blue), organ (orange) and lesion (green) levels. Circles indicate unknown quantities, i.e. model parameters to estimate, and rectangles indicate observed variables, i.e. longitudinal and survival data.

Figure S3 :

 S3 Figure S3: Time-dependent AUC (left panels) and Brier Scores (right panels) for each landmark time; s = 0 months, s = 3 months, s = 6 months, s = 9 months, and s = 12 months (from lightest to darkest respectively), depending on the joint model; the tumor burden model (solid lines), the tumor burden and range model (dotted lines) and the organ tumor burden model (dashed lines).

Table S1 :

 S1 TableS1shows the classic split-R and effective sample sizes for each candidate model. Overall, the organ tumor burden model provided the largest effective sample sizes and the R closest to 1, suggesting a better convergence for this model. The convergence was nonetheless satisfying for the tumor burden model and the tumor burden and range model. Convergence diagnostics for each candidate models; classic split R and effective sample size (neff)

		Tumor burden model Tumor burden and range model Organ tumor burden model
	Parameters	R	neff	R	neff	R	neff
	µ T S0	1.00	954	1.00	1081	1.00	759
	µ ϵ	1.00	539	1.01	403	1.00	618
	µ g	1.00	460	1.01	420	1.00	506
	µ c	1.00	595	1.00	651	1.00	669
	ξ lymph,T S0	1.00	813	1.00	1013	1.00	837
	ξ lymph,ϵ	1.01	580	1.01	458	1.00	557
	ξ lymph,g	1.00	610	1.00	584	1.00	521
	ξ lymph,c	1.00	523	1.01	695	1.00	690
	ξ lung,T S0	1.00	952	1.00	1093	1.00	649
	ξ lung,ϵ	1.00	818	1.00	654	1.00	688
	ξ lung,g	1.01	396	1.02	318	1.01	397
	ξ lung,c	1.00	319	1.02	198	1.00	287
	ξ liver,T S0	1.00	935	1.00	1290	1.00	733
	ξ liver,ϵ	1.00	731	1.00	670	1.00	728
	ξ liver,g	1.00	527	1.01	435	1.00	606
	ξ liver,c	1.00	489	1.00	727	1.00	602
	ω 1,T S0	1.00	580	1.00	543	1.01	526
	ω 1,ϵ	1.00	538	1.00	543	1.00	607
	ω 1,g	1.02	237	1.00	250	1.01	210
	ω 1,c	1.01	255	1.01	267	1.01	235
	ω 2,T S0	1.01	411	1.01	434	1.00	490
	ω 2,ϵ	1.01	552	1.00	573	1.01	612
	ω 2,g	1.00	257	1.04	125	1.00	227
	ω 2,c	1.00	206	1.02	168	1.01	203
	σ 1,lymph	1.00	521	1.00	312	1.00	496
	σ 1,lung	1.00	513	1.00	634	1.00	571
	σ 1,liver	1.00	483	1.00	486	1.00	593
	σ 1,bladder	1.00	708	1.00	506	1.00	580
	σ 2,lymph	1.00	475	1.01	270	1.00	499
	σ 2,lung	1.00	416	1.00	320	1.00	368
	σ 2,liver	1.00	308	1.00	415	1.00	485
	σ 2,bladder	1.00	433	1.03	163	1.00	455
	γ	1.00	774	1.00	672	1.00	855
	λ	1.00	612	1.00	725	1.00	806
	β	1.00	685	1.00	600	-	-
	β range	-	-	1.00	566	-	-
	β lymph	-	-	-	-	1.00	622
	β lung	-	-	-	-	1.00	579
	β liver	-	-	-	-	1.00	873
	β bladder	-	-	-	-	1.00	744

Table S2 :

 S2 Results are consistent across all candidate models. Estimates of population parameters of the tumor burden model, posterior mean [95% credibility interval] of the fixed-effect, inter-patient variability (IPV), inter-lesion variability (ILV) and location fixed-effects for each organ.

	Parameters	Baseline tumor	Treatment induced	Natural tumor	Resistance
		size T S 0 (mm)	shrinkage ϵ (day -1 )	growth g (day -1 )	appearance c (day -1 )
	Fixed-effect µ	27.0 [25.8;28.1]	0.0024 [0.0016;0.0034] 0.0031 [0.0025;0.0037] 0.0018 [0.0008;0.0031]
	IPV ω 1	0.25 [0.21;0.29]	1.31 [1.04;1.63]	0.83 [0.63;1.02]	1.36 [0.80;2.14]
	ILV ω 2	0.36 [0.33;0.39]	0.69 [0.55;0.86]	0.30 [0.06;0.48]	0.78 [0.24;1.27]
	Location effect ξ
	Lymph	-0.24 [-0.28;-0.18]	0.43 [0.18;0.70]	-0.44 [-0.64;-0.23]	1.28 [0.72;2.04]
	Lung	-0.18 [-0.24;-0.12]	0.15 [-0.15;0.44]	0.06 [-0.16;0.26]	0.41 [-0.46;1.39]
	Liver	0.02 [-0.05;0.08]	-0.31 [-0.65;0.05]	0.62 [0.42;0.82]	-1.56 [-2.46;-0.71]
	Bladder	0.40 [0.32;0.48]	-0.28 [-0.73;0.13]	-0.24 [-0.53;0.029]	-0.12 [-2.11;1.07]
	Parameters	Baseline tumor	Treatment induced	Natural tumor	Resistance
		size T S 0 (mm)	shrinkage ϵ (day -1 )	growth g (day -1 )	appearance c (day -1 )
	Fixed-effect µ	27.0 [25.7;28.3]	0.0024 [0.0015;0.0033] 0.0031 [0.0024;0.0038] 0.0018 [0.0010;0.0032]
	IPV ω 1	0.25 [0.20;0.30]	1.32 [1.06;1.64]	0.86 [0.67;1.06]	1.35 [0.77;2.15]
	ILV ω 2	0.36 [0.33;0.39]	0.67 [0.54;0.83]	0.29 [0.05;0.50]	0.76 [0.09;1.27]
	Location effect ξ
	Lymph	-0.24 [-0.30;-0.18]	0.44 [0.17;0.71]	-0.42 [-0.61;-0.23]	1.24 [0.70;1.89]
	Lung	-0.18 [-0.24;-0.11]	0.18 [-0.12;0.48]	0.06 [-0.17;0.26]	0.45 [-0.43;1.55]
	Liver	0.02 [-0.05;0.08]	-0.32 [-0.66;0.03]	0.60 [0.40;0.81]	-1.55 [-2.48;-0.70]
	Bladder	0.40 [0.32;0.48]	-0.29 [-0.78;0.12]	-0.24 [-0.50;0.035]	-0.15 [-2.01;1.06]
	Table S3: Estimates of population parameters of the tumor burden and range model, posterior mean [95% credibility
	interval] of the fixed-effect, inter-patient variability (IPV), inter-lesion variability (ILV) and location fixed-effects for
	each organ.
	Parameters	Baseline tumor	Treatment induced	Natural tumor	Resistance
		size T S 0 (mm)	shrinkage ϵ (day -1 )	growth g (day -1 )	appearance c (day -1 )
	Fixed-effect µ	27.0 [25.9;28.2]	0.0024 [0.0016;0.0035] 0.0031 [0.0025;0.0037] 0.0018 [0.0010;0.0031]
	IPV ω 1	0.25 [0.20;0.30]	1.29 [1.05;1.58]	0.82 [0.63;1.02]	1.29 [0.79;1.92]
	ILV ω 2	0.36 [0.33;0.39]	0.67 [0.53;0.81]	0.28 [0.09;0.47]	0.81 [0.26;1.24]
	Location effect ξ
	Lymph	-0.24 [-0.30;-0.18]	0.43 [0.16;0.70]	-0.41 [-0.61;-0.21]	1.22 [0.62;1.86]
	Lung	-0.18 [-0.24;-0.12]	0.20 [-0.11;0.52]	0.10 [-0.12;0.30]	0.42 [-0.35;1.37]
	Liver	0.02 [-0.05;0.09]	-0.32 [-0.64;0.013]	0.58 [0.38;0.78]	-1.54 [-2.43;-0.68]
	Bladder	0.40 [0.32;0.48]	-0.31 [-0.75;0.075]	-0.26 [-0.53;0.019]	-0.10 [-1.77;0.98]

Table S4 :

 S4 Estimates of population parameters of the organ tumor burden model, posterior mean [95% credibility interval] of the fixed-effect, inter-patient variability (IPV), inter-lesion variability (ILV) and location fixed-effects for each organ.

January 2022 of the lung lesions was much more modest (β lung 95% CI: [0.0033;0.010]).

In combination with tumor burden, the association parameter quantifying the impact of the range between the maximum and the minimum of the lesions on the instantaneous risk of death was significantly negative (β range 95% CI: [-0.013;-0.0013]). Consequently, for a given sum of individual lesions, a greater intra-patient variability induced a more favorable survival probability.

When comparing the WAIC criteria of the candidate models, combining the tumor burden impact with a marker of the intra-patient heterogeneity (tumor burden and range model WAIC=17804) or considering an heterogeneous impact of the individual lesions depending on their location (organ tumor burden model WAIC=17803) improved the fit of the model to the data as compared to the model neglecting the heterogeneity across lesions (tumor burden model WAIC=17810). Alternative association structures were considered, such as the maximum of the lesions, but they provided a poorer adjustment to the data (see Table S5 in supplementary material).

Goodness-of-fit

All candidate models provided satisfying population parameters estimates: convergence was systematically achieved and all population parameters R were below 1.04 overall, and below 1.01 in the organ tumor burden model. The median (minimum-maximum) effective sample size over the population parameters was 527 in the tumor burden model, 524 in the tumor burden and range model, and 586 (203-873) in the organ tumor burden model. In summary, the tumor burden and range model had slightly deteriorated estimation precision and stability of the population parameters, especially regarding the inter-lesions variability of the biological parameters g and c (see Table S1 in Supplementary material).

T., Hasegawa, T., Uchibori, K., Yanagitani, N., Horiike, A., Horai, T., Seike, M., Gemma, A., and Nishio, M. (2020). Dissociated responses at initial computed tomography evaluation is a good prognostic factor in non-small cell lung cancer patients treated with anti-programmed cell death-1/ligand 1 inhibitors. BMC Cancer 20, 207. Zecchin, C., Gueorguieva, I., Enas, N. H., and Friberg, L. E. (2016). Models for change in tumour size, appearance of new lesions and survival probability in patients with advanced epithelial ovarian cancer. British Journal of Clinical Pharmacology 82, 717-727. Number: 3. 3 Alternative association structures

Table S5 shows the link parameter posterior mean and its 95% CI, as well as the WAIC criterion for the alternative candidate models. As expected, the maximum of the lesions is significantly associated with the risk of death when considered alone. Similarly, the range of the lesions was significantly associated with the risk of death. However, the WAIC of these models were much larger than the WAIC of the tumor burden model, the tumor burden and range model or the organ tumor model (WAIC=17810, 17804 and 17803 respectively). These results suggest the necessity of considering all the lesions together, instead of just the maximal lesion or the range alone.

Model

Definition of the reference joint model of SLD and survival

We compared the predictive performances of our final individual lesions model with a reference model of the Sum of the Longest Diameters of the target lesions and survival. In this reference model, only the measurements of the SLD of the target lesions were fitted for the longitudinal part. Denoting z i,l the l th of SLD for patient i measured at time t i,l , and T S the sTGI structural model defined in section 2.2:

Formally, for each biological parameter q ∈ {T S 0 , ϵ, g, c}, individual parameter ψ i,q is function of a population fixed-effect µ q and a patient random effect normally distributed η i,q ∼ N (0, ω 2 1,q ): log(ψ i,q ) = log(µ q ) + η i,q .

We assumed an association between the current value of the SLD and the instantaneous risk of death:

with h 0 (t) = γ λ ( t λ ) γ-1 a Weibull baseline hazard function with shape parameter γ and scale parameter λ, T S(t, ψ i ) the SLD at time t predicted by the model, and β SLD the link parameter characterizing the association between SLD and survival.

Estimates of the parameters of the reference joint model of SLD and survival

We reported in table S6 the estimates of longitudinal parameters of the SLD and survival joint model. In this model, we do not take into account the impact of tumor location neither the ILV. Hence, the IPV is inflated as compared to the IPV in the joint model of individual lesions and survival, especially for the tumor growth (g) and the resistance appearance (c) parameters.

Fixed-effects

Time-dependent Brier Scores

In Figure S2, we shows the time-dependent Brier Scores under the organ tumor burden model and the SLD model for each landmark times. At landmark s = 0 and s = 90 days, the organ tumor burden model showed a better discrimination ability in comparison to the SLD model in the total population. This phenomenon was exacerbated in the subpopulation of patients having liver or bladder lesions.

Comparison of predictive performances according to the association structure

In figure S3, we compared time-dependent AUC and Brier Scores for the three main association structures, namely the tumor burden model, the tumor burden and range model, and the organ tumor burden model. Model predictive


Benjamin Grain

Benoît Snella

Jérôme Vicogne

Frédéric Capet

Birgit Wiltschi

Oleg Melnyk

Vangelis Agouridas

Fast-SEA ligation: a chemoselective NCL-inspired ligation reaction effective in the nanomolar concentration range

Introduction

The Chemical Protein Synthesis has highly benefitted from the properties of the thioester functional group, in particular with the development of the Native Chemical Ligation (NCL) reaction. 1,2 However, NCL-based modification of proteins suffers from the moderate reactivity of classical alkyl and aryl thioesters 3 , especially in highly dilute media. Seeking for fast reactive thioesters surrogates, we have discovered that oxalyl amide SEA ( oxo SEA), in the form of a latent thioester precursor, 4 could, upon activation, act as a powerful acylating agent in chemoselective ligation reactions with β-aminothiols such as the cysteine residue. Investigation on the reaction of oxo SEA peptide 2 with cysteinyl peptide 3.

To further explore the performance of the oxo SEA-mediated ligation in the nanomolar concentration range, we opted for SDS-PAGE separation followed by Western Blot as the analysis tool. This allowed us to work at the protein level by monitoring the reaction between the oxo SEA peptide derived from V5 peptidic tag 5 and the C-terminally biotinylated kringle 1 domain of the hepatocyte growth factor (K1/HGF) equipped with a N-terminal cysteine 6. The groundbreaking reactivity of this activatable oxalyl thioester surrogate allows us to perform peptide and protein modification through a chemoselective ligation reaction. The latter proceeds extremely fast in aqueous media, under mild conditions and in highly dilute medium which can be as complex as a cell lysate. By being easy to introduce within peptide segments in the form of a lysine derivative, oxalyl thioesters might be of significant importance for the development of novel crosslinking strategies or the polycondensation of water-soluble (bio)molecules. Moreover, this ligation reaction is also efficient and chemoselective in a complex protein medium, as it has been successfully performed in a crude protein extract of a THP-1 whole cell lysate supplemented with 100 nM of 6 in the presence of a 15-fold excess of 5 (see article).

X-ray structure of 1

( 1 )

 1 Dawson, P. E. et al, Science 1994, 266, 776. (2) Agouridas, V. et al, Chem. Rev. 2019, 12, 7328. (3) Saito, F. et al, ACS Chem. Biol. 2015, 10, 1026.

  the gram scale with 17% overall yield -Introduction of Fmoc-Lys( oxo SEA)-OH 1 within peptide segments at any position, according to standard Fmoc-SPPS protocols a Univ. Lille, CNRS, Inserm, CHU Lille, Institut Pasteur de Lille, U1019 -UMR 9017 -CIIL -Center for Infection and Immunity of Lille, F-59000 Lille, France b Institute of Bioprocess Science and Engineering, University of Natural Resources and Life Sciences, Vienna c Centrale Lille, F-59000, Lille, France

-

  a) Ollivier, N. et al, Org. Lett. 2010, 12, 5238. b) Diemer, V. et al, in Total Chemical Synthesis of Proteins (Eds.: A. Brik, L. Liu, P. Dawson), Wiley, 2021, 59.oxo SEA ligation with a peptide concentration of 1 mM is: -Efficient (HPLC conversions to polypeptide 4 beyond 95%) -Fast (Half-reaction time of 200 s) -Optimal at a working pH = 5,5 -Triggered by reductants (TCEP) rate constant -Rate break at 150 µM, and below this concentration value, apparent second-order rate constant thrives around 30 M -1 s -1 Detection of ligated product 7 after 200 minutes of reaction oxo SEA ligation effective in highly dilute solution, with peptide concentration of 270 nM

Acknowledgements to the Agence Nationale de la Recherche (ANR-21-CE44-0031) and the University of Lille (R-PILOTE-19-0008-Molecular).

More details on this study: Snella, B. et al, ChemRxiv.

2022.

Comparison with classical NCL rates highlighted the tremendous reactivity of the oxo SEA group, even more when decreasing peptide concentration to 50 µM.


Benoît Snella

Benjamin Grain

Jérôme Vicogne

Frédéric Capet

Birgit Wiltschi

Oleg Melnyk
email: oleg.melnyk@ibl.cnrs.fr

Vangelis Agouridas
email: vangelis.agouridas@ibl.cnrs.fr

Fast protein modification in the nanomolar concentration range using an oxalyl amide as latent thioester

We show that latent oxalyl thioester surrogates are a powerful means to modify peptides and proteins in highly dilute conditions in purified aqueous media or in mixtures as complex as cell lysates. Designed to be shelf-stable reagents, they can be activated on-demand for enabling ligation reactions down to peptide concentrations as low as a few hundreds nM at rates approaching 30 M -1 s -1 .

Thioester based-chemistry has been the source of many innovative breakthroughs in domains such as chemical biology 1 or materials science. [START_REF] Aksakal | Thioester functional polymers[END_REF][START_REF] Worrell | A user's guide to the thiol-thioester exchange in organic media: scope, limitations, and applications in material science[END_REF] In particular, the reaction of a thioester group located at the C-terminus of peptides or proteins with a cysteinyl (Cys) peptide is the cornerstone of the Native Chemical Ligation (NCL) reaction [START_REF] Dawson | Synthesis of proteins by native chemical ligation[END_REF] and related methods (Figure 1a). [START_REF] Agouridas | Native Chemical Ligation and Extended Methods: Mechanisms, Catalysis, Scope, and Limitations[END_REF] The hallmark of this process unquestionably resides in its exquisite chemoselectivity and its capacity to take place in water, under very mild pH and temperature conditions. During the last three decades, this reaction has profoundly remodeled protein chemistry by enabling the routine production of homogeneous batches of native or modified (semi-)synthetic proteins [START_REF] Conibear | Native chemical ligation in protein synthesis and semi-synthesis[END_REF][START_REF] Muir | Expressed protein ligation: a general method for protein engineering[END_REF] . [START_REF] Agouridas | A statistical view of protein chemical synthesis using NCL and extended methodologies[END_REF] Considering the ease with which the Cys residue can be introduced at the N-terminus of peptides and proteins through chemical approaches or recombinant techniques, NCL is also an appealing reaction for accessing protein conjugates. [START_REF] Boutureira | Advances in Chemical Protein Modification[END_REF] For the latter application, achieving rapid reactions with N-terminal Cys is highly desirable to accommodate protein limited solubility and minimize side-reactions. Several kinetic studies have shown that the archetypal NCL reaction performed with alkyl or aryl thioesters generally exhibits slow to moderate rates under standard experimental conditions (~0.3 to 4 M -1 s -1 ). [START_REF] Diemer | A cysteine selenosulfide redox switch for protein chemical synthesis[END_REF][START_REF] Saito | Critical Evaluation and Rate Constants of Chemoselective Ligation Reactions for Stoichiometric Conjugations in Water[END_REF] While powerful in the mM concentration range, the NCL reaction shows its limits at µM concentrations due to competitive thioester hydrolysis. [START_REF] Bracher | The relative rates of thiol-thioester exchange and hydrolysis for alkyl and aryl thioalkanoates in water[END_REF] Thus, the use of NCL-based approaches becomes challenging for any application proceeding at high dilution such as the modification of proteins, although selenoester chemistry offers some opportunities but to the expense of increased complexity. [START_REF] Chisholm | Peptide Ligation at High Dilution via Reductive Diselenide-Selenoester Ligation[END_REF] In search for fast reacting thioester derivatives, our attention was attracted by oxalyl thioesters. These species were first reported as highly reactive and transient mammalian metabolites, which could be trapped in cell lysate with cysteine following an NCL-like mechanism (Figure 1b). [START_REF] Skorczynski | Oxalyl thiolesters and N-oxalylcysteine are normal mammalian metabolites[END_REF] Their reactivity in water remains however largely unexplored, probably due to the difficulties encountered in installing them on biomolecules and to their ease of hydrolysis with a half-lives of a few minutes in water at neutral pH. 15 Against the odds, we discovered that oxalyl amides in the form of bis(2-sulfanylethyl)amido (SEA [START_REF] Ollivier | Bis(2sulfanylethyl)amino native peptide ligation[END_REF][START_REF] Diemer | Se-acyl transfer devices in protein synthesis[END_REF] ) latent thioester precursors, enable an NCL-like reaction to proceed with exceptional rates approaching 30 M -1 s -1 (Figure 1c). This result is unexpected because a feature of such rearranging amides is the modest rate of the N,S-acyl group migration step, making them unsuitable as efficient acyl-transfer agents in highly dilute conditions. [START_REF] Agouridas | Native Chemical Ligation and Extended Methods: Mechanisms, Catalysis, Scope, and Limitations[END_REF] We therefore undertook an in-depth investigation of the reactivity of the oxo SEA system, which enabled protein modification to proceed efficiently in the nanomolar concentration range in purified reaction medium or in mixtures as complex as cell lysates.

The model system used for investigating the reactivity of SEA-based oxalyl thioester precursors is presented in Figure 2a. While the incorporation of an oxo SEA group at the N-terminus of peptides during Solid Phase Peptide Synthesis is a straightforward method for rapid functionalization (see Supplementary Information), we have opted for a more flexible approach by preparing a conveniently protected lysine derivative 1, i. e. Fmoc-Lys( oxo SEA)-OH, which has the advantage of being easily introduced at any position into peptide sequences. Lysine derivative 1 was obtained at the gram scale after 6 steps from commercially available starting materials and with 17% overall yield. Its structure was validated by X-ray diffraction crystallography (Figure 2b). Amino acid 1 was easily incorporated into peptides of type 2 using standard Fmoc-SPPS procedures and proved perfectly stable during peptide cleavage and purification steps. In a first series of experiments, we used model peptides 3a and 4a to evaluate the potential of oxo SEA-containing peptides as acyl-transfer agents in different experimental conditions using MPAA [START_REF] Johnson | Insights into the mechanism and catalysis of the native chemical ligation reaction[END_REF] as the thiol catalyst [START_REF] Diemer | Pedal to the Metal: The Homogeneous Catalysis of the Native Chemical Ligation Reaction[END_REF] and TCEP as the reductant (Figure 2c). The presence of a tryptophan residue in peptide 3a enabled accurate quantification of all the species of interest on HPLC chromatograms by UV-monitoring at 280 nm, including at the lowest peptide concentrations tested. The corresponding data were used for the determination of rate constants by nonlinear regression fitting.

Preliminary experiments established the influence of additives concentration and of the pH on ligation efficiency at 1 mM peptide concentration (see Supplementary Information). The reaction was found to be dependent on the concentration of MPAA but still displayed acceptable conversion ratios in the absence of a thiol catalyst. Importantly, the oxo SEA group remained stable and silent in the absence of TCEP or at low concentrations of MPAA, giving the possibility to trigger it on-demand. These experiments revealed also that the oxo SEA group is reacting significantly faster at the mildly acidic pH of 5.5 than at neutral pH. In that, the oxo SEA group apparently follows the behavior of C-terminal SEA peptides for which mildly acidic conditions are known to promote the N,S-acyl shift process. [START_REF] Cargoët | Catalysis of thiol-thioester exchange by water-soluble alkyldiselenols applied to the synthesis of peptide thioesters and SEA-mediated ligation[END_REF] Therefore, we initially performed ligations at pH 5.5. The reaction of peptides 3a and 4a, respectively introduced at 1 mM and 1.2 mM, yielded the expected ligated product 5a with an astonishing half-reaction time (t1/2) of 200 s for a final HPLC yield >95% (Figure 3a). For comparison, classical SEA ligations performed in the same concentration range require about 5-7 hours to be half-completed, highlighting the extreme difference in reactivity of the SEA and oxo SEA groups. Gradually, lowering the concentration of peptides 3a and 4a to 50 µM revealed that the reaction was still proceeding very efficiently (t1/2 = 900 s). Ligation product 5a was repeatedly obtained with a conversion ratio close to 90%, showing that even used at low molarities, the highly reactive oxalyl thioester surrogates were not significantly hydrolyzed. As a control, NCL reactions performed with classical alkyl and aryl thioesters at the same peptide concentration proceeded more slowly and plateaued at ~65-70% conversion due to hydrolysis of the starting material (Figure 3b). At 5 µM of peptide 3a and still at pH 5.5, the hydrolysis of the oxalyl thioester became noticeable, thus reducing the yield of ligated peptide 5a (~40%). At this point, we reasoned that at such a low concentration the rearrangement of the amide might not be limiting anymore. If this is the case, the optimal pH conditions identified from reactivity studies conducted at mM peptide concentration might no longer apply to dilute conditions. In particular, working at neutral pH would perhaps be more appropriate for promoting the thiolthioester exchanges and thus the capture of the oxalyl thioester by the Cys peptide. [START_REF] Agouridas | Native Chemical Ligation and Extended Methods: Mechanisms, Catalysis, Scope, and Limitations[END_REF] This proved to be the case since a slight pH adjustment from 5.5 to 7.2 was sufficient to restore satisfying levels of conversion at 5 µM (80%, Figure 3a). Intriguingly, the determination of apparent second order rate constants for the oxo SEA-mediated process at various peptide concentrations revealed an uncommon behavior (Figure 3c). While the rate constant obtained at 1 mM was ~4 M -1 s -1 , it counterintuitively increased as peptide concentration decreased, with a remarkable break at 150 µM. Below this concentration threshold, the apparent second order rate constant stabilized at a high level of ~30 M -1 s -1 .

UPLC-MS monitoring of the reaction mixture at various peptide concentrations furnished useful information regarding the nature of some reaction intermediates involved in this process (Figure 3d).The observation of reduced oxo SEA species 6 and rearranged thioester 7 are expected according to previous studies on the C-terminal SEA group in peptides (Figure 3e). [START_REF] Cargoët | Catalysis of thiol-thioester exchange by water-soluble alkyldiselenols applied to the synthesis of peptide thioesters and SEA-mediated ligation[END_REF] The oxo SEA-derived thioester 7 likely undergoes a thiol-thioester exchange with the aryl thiol catalyst MPAA to form 8, which could not be observed in our experiments for being probably quickly consumed by the Cys peptide to produce the desired ligated product 5.

Here, the peculiarity arises from the formation of an unusual dimeric intermediate 9 observed in UPLC. This adduct significantly accumulates within the mixture when the reaction is conducted at 1 mM peptide concentration but falls to basal levels at 100 µM and below. The accumulation of dimer 9 at high peptide concentration shows that it is slowly reacting. It must first dissociate to evolve into ligated product 5 and its formation likely accounts for the dependence of the rate constant to peptide concentration. Importantly, ligation proceeds extremely fast at low peptide concentration, that are the conditions for which the chemical system is primarily intended for.

Thus far, the results obtained for the oxo SEA-mediated ligation on model peptides 3a and 4a provided suggestive evidence that the reaction might still operate in the nanomolar concentration range. However, 5 µM is the lowest working concentration limit allowing reliable analytical detection and reproducible determination of ligation yields by HPLC. To further explore the performance of the reaction in the nanomolar concentration range, we opted for SDS-PAGE separation followed by Western-Blot as the analysis tool and to work at the protein level. For this purpose, we produced the oxo SEA peptide derived from V5 peptidic tag 3b and the 10 kDa C-terminally biotinylated kringle 1 (K1) domain of the hepatocyte growth factor (K1/HGF) equipped with an N-terminal Cys 4b. [START_REF] Ollivier | A simple and traceless solid phase method simplifies the assembly of large peptides and the access to challenging proteins[END_REF] V5 tag can be detected with high sensitivity using a specific anti-V5 monoclonal antibody labelled with horseradish peroxidase (HRP), while streptavidin-HRP enables to trace the K1 protein domain. The latter was reacted at concentrations ranging from 2 µM to 270 nM in the presence of ~6 equivalents of V5 tag oxo SEA peptide 3b. The analysis shown in Figure 4 allows to appreciate the successful formation of the ligation product 5b after 200 minutes, even when K1 polypeptide 4b was introduced at a concentration of 270 nM. Importantly, the control experiments performed in the absence of TCEP or in the presence of a biotinylated K1 domain lacking the N-terminal cysteine residue failed to produce any ligation product, respectively demonstrating the latency of the system and its chemoselectivity (see Supplementary Information). After having characterized the oxo SEAmediated ligation in a purified media, we were interested to know if the same could operate in mixtures as complex as a cell lysate (see Supplementary Information). To this end, the protein extract of a THP-1 whole cell lysate was supplemented in Cys-K1B 4b at 100 nM final concentration and reacted in the presence of 15 equivalents of V5-Lys( oxo SEA) peptide 3c. Strikingly, a band revealed by both anti-V5 antibody and streptavidin-HRP conjugate at ~15 kDa suggested the formation of the ligation product 5c. Its formation was however accompanied by an important background noise and the appearance of other bands. Controls at different reaction times, in the presence or in the absence of a strong reducing agent enabled to conclude that a nonspecific interaction of peptide 3c introduced in excess with unidentified elements of the protein extract accounted for this observation. The successful formation of ligated product 5c in such a complex mixture again highlighted the high selectivity and efficiency of the oxo SEA-mediated ligation for protein modification.

In conclusion, we herein report the exceptional reactivity of a bioinspired and activatable oxalyl thioester surrogate for peptide and protein modification at high dilution. This highly chemoselective reaction proceeds remarkably fast in aqueous media, under mild conditions and allows to perform ligations in purified and highly complex media such as cell lysates. Beyond protein modification, this work extends the scope of oxalyl thioesters as valuable reagents by addressing the issue of their poor stability in water through the synthesis of bench-stable and activatable surrogates. In doing so, these compounds or any derivative thereof might be used for future works as a complementary approach for practical and widely applicable crosslinking strategies or polycondensation of water-soluble (bio)molecules.

Figure 1 .

 1 Figure 1. a) General principle of the NCL reaction. b) Oxalyl thioester metabolites are prone to

Figure 2 .

 2 Figure 2. a) Synthesis of the Fmoc-Lys( oxo SEA)-OH residue 1 and its incorporation into peptides of type 2. b) ORTEP diagram of the X-ray structure of Fmoc-Lys( oxo SEA)-OH 1 with ellipsoids shown at the 50% probability level. c) Model oxo SEA-mediated ligations examined in this study.

Figure 3 .

 3 Figure 3. a) Conversion to ligated peptide 5a from peptides 3a and 4a. Unless otherwise specified,

Figure 4 .

 4 Figure 4. a) oxo SEA-mediated ligation between V5 containing peptide 3b or 3c and biotinylated K1

Acknowledgements

This research was supported by the Agence Nationale de la Recherche (ANR-21-CE44-0031), the University of Lille (R-PILOTE-19-0008-Molecular) and the Austrian Science Fund (FWF I 5800-N). We thank CNRS, INSERM, Institut Pasteur de Lille and University of Lille for core facilities. Hervé Drobecq is gratefully acknowledged for the proteomic analysis of the ligation product 5a. Keywords: non-native ligation • thioester • protein modification • N,S-acyl shift systems • latency