Bayesian modelling of individual lesions dynamics and survival to characterize response to immunotherapy cancer treatments

Kerioui M (1), (2), (3), (4), Bertrand J (1), Mercier F (5), Bruno R (4), (6), Desmée S (2), Guedj J (1)

(1) Université de Paris, INSERM, IAME, F-75018 Paris, France (2) Université de Tours, Université de Nantes, INSERM SPHERE, UMR 1246, Tours, France (3) Clinical Pharmacology, Roche/Genentech, Paris, France (4) Institut Roche, Boulogne-Billancourt, France (5) Biostatistics – Roche Innovation Center Basel, Basel Switzerland (6) Clinical Pharmacology, Genentech Inc., Marseille, France

Context: In oncology clinical trials, a common marker of tumor size is the Sum of Longest Diameters (SLD) of target lesions, and is part of the RECIST criteria, which evaluate response and progression. To characterize the association between the tumor size dynamics and survival, one can simultaneously fit both longitudinal and time-to-event data using a joint model.

With this approach, early information on response to treatment contained in patient tumor dynamics might help to identify most at risk patients[1], thus, optimize clinical trial[2] and improve patients follow-up. However, using a composite marker such as SLD has several limitations, as it neglects the heterogeneity in lesion dynamics, which might be partly explained by their location[3]. This might be even more damaging under immunotherapy treatments, as they might increase the risk of dissociated responses across lesions, especially in different organs[4]. As immunotherapies are currently revolutionizing the treatment of cancer patients, there is an urgent need to understand the exacerbated variability in response to those treatments[5].

Incorporating all the source of variability into one model might help improving its prediction abilities, but raises many methodological questions; how to differentiate tumor kinetics and its impact on survival in each organ, how to characterize the impact of responding and progressing lesions within one same patient? For such a hierarchical structure and the computational challenge it involves, Bayesian inference could be more efficient. In practice, the Hamiltonian Monte-Carlo (HMC) algorithm is known for its good inference properties for complex models[6], but has never been assessed for nonlinear joint model inference.

To address these questions, we relied on a rich dataset from a phase 2 (IMvigor210[7]) and a phase 3 (IMvigor211[8]) clinical trials of 300 and 900 advanced bladder cancer patients, treated with atezolizumab. Each patient had up to 5 target lesions followed over time, resulting in a total of 2133 target lesions and 6697 measurements in the phase 3.

Objectives:

1. To evaluate the HMC algorithm implemented in Stan software for nonlinear joint model parameters Bayesian inference.
2. To develop a model integrating the effects of tumor location on response dynamics and survival and evaluate its benefit to predict a clinical trial outcome.
3. To develop a multilevel joint model to quantify the inter-lesions variability under immunotherapy and investigate its benefit to improve individual survival prediction.

Methods:

1. We aimed to validate an innovative Bayesian algorithm for nonlinear joint model parameters inference. HMC inference properties were first assessed in a simulation study. Then, we analyzed IMvigor210 data under three different prior scenarios. Using a cross-validation method[9], we selected the best association structure between SLD kinetics and survival among four candidates. We assessed the goodness-of-fit of the model based on individual fits, and posterior predictive checks.
2. We built a model integrating the effects of organ-specific tumor size on survival in IMvigor211 data, considering the five following locations: lymph, lung, liver, bladder and other sites. We assessed the benefit of our model in predicting the trial outcome based on early data compared with a SLD-based model. We re-estimated model parameters, assuming early cutoffs of tumor size and survival data after last patient inclusion. Then, we emulated 1000 trials and we evaluated trial outcome (i.e. atezolizumab over chemotherapy hazard ratio) in each replicate.
3. We developed a multilevel joint model to describe individual target lesions and their impact on survival in patients from IMvigor211 immunotherapy arm. We used a hierarchical structure including a first usual patient level, a second organ level, and a third lesion level. Dynamic prediction approaches were used to assess the benefit of individual lesions follow-up in detecting most at risk patients, compared to total SLD follow-up.

Results:

1. HMC provided unbiased estimates of population parameters and precisely handled uncertainty with satisfying 95% coverage rates. Interestingly, the link parameter, that characterizes the association between tumor size and survival, was precisely estimated, with a coverage rate close to 95%. Results were confirmed in real data analysis, and our model showed good adjustment to the data regardless of prior information[10]. Of note, Bayesian inference was computationally intensive compared to frequentist approach, and the latter is therefore preferred for simple models.
2. The Stochastic Approximation of Expectation-Maximization algorithm in Monolix software performed well for population parameters estimation. However, we approached the limits of frequentist inference, with large standard errors on some parameters. We found significantly different tumor dynamics from one organ to another, notably a faster growth of the liver tumor, and a stronger treatment effect in the lymph nodes. We showed the importance of integrating the impact of tumor location on the prediction of survival. As an illustration, patients having a 10 mm growth of their liver lesion had their instantaneous risk of death increased by 12%. In contrast, the same growth in the lung increased the risk of death by only 5%. The treatment effect significance could be predicted based on a 6-months tumor size and survival cutoff with the organ-specific tumor size joint model when a longer follow-up (9 months) was required with the usual SLD-based joint model.
3. We developed three multilevel joint model in Stan, describing the three different tumor size markers (total SLD, organ SLD and individual target lesions) and their association with survival, integrating up to three levels of random effects. The 95% credibility intervals of inter-lesions variability parameters excluded zero, suggesting a non-negligible within-patient variability. About 16% of patients treated with atezolizumab had dissociated responses, defined as a 30%-decrease from baseline and a 20%-increase above nadir in two different lesions. Preliminary results suggested that for these patients, the lesion-level model was particularly adapted to provide a more precise survival prediction, and should be confirmed by ongoing dynamic prediction analysis.

Conclusion: We first showed that HMC algorithm could provide precise estimation of nonlinear joint model parameters. Second, we showed the benefit of organ-specific tumor size follow-up in understanding the underlying mechanism of the response to immunotherapy, and predicting the outcome of a trial based on early data. Finally, we developed a joint model to describe individual lesions and their association with survival. We are currently comparing the benefit of this model over less detailed model, to improve the early detection of patients in need of specific follow-up under immunotherapy.

References:
[1] Tardivon C, Desmée S, Kerioui M, et al. Association Between Tumor Size Kinetics and Survival in Patients With Urothelial Carcinoma Treated With Atezolizumab: Implication for Patient Follow-Up. Clin Pharmacol Ther. 2019;106(4):810-820. doi:10.1002/cpt.1450
[2] Claret L, Jin JY, Ferté C, 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. 2018;24(14):3292-3298. doi:10.1158/1078-0432.CCR-17-3662
[3] Mercier F, Kerioui M, Desmée S, Guedj J, Krieter O, Bruno R. Longitudinal analysis of organ-specific tumor lesion sizes in metastatic colorectal cancer patients receiving first line standard chemotherapy in combination with anti-angiogenic treatment. J Pharmacokinet Pharmacodyn. Published online 31 August 2020. doi:10.1007/s10928-020-09714-z
[4] Humbert O, Chardin D. Dissociated Response in Metastatic Cancer: An Atypical Pattern Brought Into the Spotlight With Immunotherapy. Front Oncol. 2020;10:566297. doi:10.3389/fonc.2020.566297
[5] Borcoman E, Kanjanapan Y, Champiat S, et al. Novel patterns of response under immunotherapy. Ann Oncol. 2019;30(3):385-396. doi:10.1093/annonc/mdz003
[6] Monnahan CC, Thorson JT, Branch TA. Faster estimation of Bayesian models in ecology using Hamiltonian Monte Carlo. O’Hara RB, ed. Methods Ecol Evol. 2017;8(3):339-348. doi:10.1111/2041-210X.12681
[7] Rosenberg JE, Hoffman-Censits J, Powles T, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. The Lancet. 2016;387(10031):1909-1920. doi:10.1016/S0140-6736(16)00561-4
[8] Powles T, Durán I, van der Heijden MS, et al. Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial. The Lancet. 2018;391(10122):748-757. doi:10.1016/S0140-6736(17)33297-X
[9] Vehtari A, Lampinen J. Bayesian Model Assessment and Comparison Using Cross-Validation Predictive Densities. Neural Computation. 2002;14(10):2439-2468. doi:10.1162/08997660260293292
[10] Kerioui M, Mercier F, Bertrand J, et al. Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy. Stat Med. Published online 8 October 2020:sim.8756. doi:10.1002/sim.8756

Reference: PAGE 29 (2021) Abstr 9628 [www.page-meeting.org/?abstract=9628]

Oral: Lewis Sheiner Student Session