Mechanistic Modeling of Joint Circulating Cell-free DNA Concentration—Tumor Size Kinetics under Immune-Checkpoint Inhibitors in Advanced Cancer
Linh Nguyen Phuong1,2,3, PhD Frédéric Fina4, MD Jean-Laurent Deville5, MD, PhD Pascale Tomasini3,5,6, MD, PhD Laurent Greillier1,2,3,5,6, MD, PhD Caroline Gaudy-Marqueste3,5,7, PhD Audrey Boutonnet8, PhD Frédéric Ginot8, PhD Jean-Charles Garcia8, MD, PhD Sébastien Salas1,2,3,5, MD, PhD Sébastien Benzekry1,2,3
1COMPO (COMPutational pharmacology and clinical Oncology), Centre Inria d'Université Côte d'Azur, 2Centre de Recherche en Cancérologie de Marseille (CRCM), Inserm U1068, CNRS UMR7258, Institut Paoli-Calmettes, 3Aix-Marseille University, 4ID-solutions oncology, 5Assistance Publique—Hôpitaux de Marseille, 6Multidisciplinary Oncology & Therapeutic Innovations Department, CEPCM, 7Dermatology and Skin Cancer Department, 8Adelis technologies
Objectives In cancer patients, the plasmatic concentration of global cfDNA — shed by both tumor and wild-type (especially immune) cells — rises and fluctuates with tumor size, number of metastases, and treatment¹. Analyzing cfDNA concentration and fragment sizes² (part of the fragmentome) is the topic of intense current research and part of the emerging field of liquid biopsy. Given the unclear biology underlying cfDNA release, fragmentation and elimination, mechanistic modeling of joint tumor-cfDNA longitudinal data can shed light on the relationship between tumor kinetics (TK), cfDNA levels, and treatment outcome. The ongoing SChISM (Size cfDNA Immunotherapies Signature Monitoring) clinical study monitors plasma cfDNA size profiles during immune-checkpoint inhibition (ICI). It aims to better understand cfDNA biology during treatment and enable early therapeutic adjustments to mitigate ICI-related progression or toxicity. The objectives of the present work were to: 1) Develop a mechanistic model of cfDNA–TK from clinical data of cancer patients undergoing ICI. 2) Study the association of early, on-treatment, model-based parameters with clinical outcome. Methods Data One hundred and thirty-nine advanced cancer patients treated with ICI, either as monotherapy or in combination therapy, were enrolled. Cancer types comprised non-small cell lung cancer (NSCLC, N = 60), squamous cell carcinomas of the head and neck (N = 28), metastatic clear cell renal cancer (N = 13), advanced urothelial bladder carcinomas (N = 9), and melanoma (N = 29). Using the patented cost-effective BIABooster³,4 device from Adelis Technologies (France), cfDNA size profiles were generated from plasma samples collected at baseline and at each drug administration. The primary endpoint was early progression (EP), defined as the progression at first imaging, around 3 months after treatment initiation. The secondary endpoint was progression-free survival (PFS). Model A mechanistic model was developed to jointly describe the dynamics of global cfDNA concentration and tumor lesions under immunotherapy. The model relied on the following biological assumptions: - Tumor cells (assimilated to the sum of largest diameters S(t)) comprise two sub-populations: cells sensitive to treatment Ss and resistant ones Sr, governed by first ordered kinetics driven by a shrinkage parameter KS and a regrowth parameter KG. - Cancer cells induce the release of DNA fragments (concentration C(t)) due to their chaotic apoptosis-proliferation cycle characteristics5, with a rate ? - CfDNA is eliminated from the circulation by a process depending only on its concentration f(C). dSs/dt = KS.Ss dSr/dt = KG.Sr S = Sr + Ss dC/dt = ?.S - f(C) with initial conditions Ss(t=0)=Ss0, Sr(t=0)=Sr0, C(t=0)=C0, where t=0 denotes treatment initiation. We compared three hypotheses for cfDNA clearance by the liver and kidneys6: - f1(C) = kD the constant clearance. - f2(C) = kD.C the proportional clearance. - f3(C) = kD.C / (C50+C) the saturated clearance. Mixed-effects modeling To account for inter-individual variability, random effects were assumed to have log-normal distribution on the individual parameters. The best error models were constant for the TK estimation and proportional for cfDNA one, which minimized the corrected Bayesian information criterion. The population parameters were identified using the Monolix software and the Stochastic-Approximation Expectation Maximization algorithm7. Tumor kinetics (TK) population parameters (Ss0, Sr0), KS, and KG) were identified from radiological report data and then fixed in the joint cfDNA-TK model. In a second step for association with PFS, the individual empirical Bayes estimates were re-identified using longitudinal data truncated at 1.5 months, using 118 patients who had not progressed. The population parameters prior was fixed from the full kinetics (FK) model for Bayesian estimation. Association with clinical outcome The truncated empirical Bayes parameter estimates were pooled with pre-treatment fragmentome-derived metrics (e.g., concentration within multiple size ranges, position of the distribution peaks) to assess associations with the endpoints, using univariable and multivariable logistic/Cox proportional hazard regression models. Results The double exponential model with constant error provided the best fit for TK. The final TK-cfDNA model assumed proportional clearance f2(C) and treated monotherapy and combination therapy equivalently. It was able to accurately describe non-trivial cfDNA kinetics. Notably, it captured both steady trends in cfDNA concentration and transient spikes at treatment initiation. These observations were impossible to generate under models with constant or saturated clearance, indicating that proportional cfDNA clearance is a likely mechanism. Goodness-of-fit diagnostics indicated no model misspecification: individual conditionally weighted residuals were centered around zero, all population parameter relative standard errors were below 19%, and parameter estimates absolute correlations did not exceed 0.22. A condition number of 3.13 further confirmed the model was not over-parameterized. In statistical FK analyses, higher KG was significantly associated with EP and shorter PFS, as well as Sr0 Both outperformed the best baseline cfDNA metric, i.e., the relative concentration of long fragments (= 1650 base pairs). CfDNA FK parameters (C0 and kD) were significantly associated with PFS only for NSCLC patients (log-rank test: p-value = 0.007 and 0.03, respectively). These results remained consistent for KG in the truncated analysis (OR UV: 1.66 [CI: 1.14-2.42], p = 0.008; AUC: 0.74; HR UV: 1.49 [CI: 1.19-1.86], p < 0.001; C-index: 0.67). KG also outperformed the best baseline cfDNA metrics, i.e., long fragments (= 1650 base pairs), conversely to S_(R_0 ). KS was also associated with longer PFS (HR: 0.31 [CI: 0.15-0.64], p = 0.003; C-index: 0.8), as well as the proportion of resistant cells at baseline (HR: 1.37 [CI: 1.05-1.77], p = 0.02; C-index: 0.63). Conclusions Our mechanistic modeling approach supports a proportional clearance of cfDNA and successfully captured cfDNA “bumps” at treatment initiation. Consistent with previous findings, early tumor kinetics parameters (growth and decay rates, as well as the baseline proportion of resistant cells) were associated with outcome. The parameters specific to cfDNA kinetics correlated with PFS only in the NSCLC subgroup and when considering FK. Future work will explore fragment size-dependent release and fragmentation mechanisms, to account for different biological mechanisms of cfDNA shedding (e.g., apoptosis for short fragments and necrosis for longer fragments6). Integrating joint modeling with survival outcomes is also an interesting avenue, alongside with multivariable analyses leveraging machine learning techniques.