Célestin Bigarré
email: celestin.bigarre@inria.fr

François Bertucci

Pascal Finetti

Gaëtan Macgrogan

Xavier Muracciole

Sébastien Benzekry

Mechanistic modeling of metastatic relapse in early breast cancer to investigate the biological impact of prognostic biomarkers

Introduction

Breast cancer is the most common cancer amongst women and has a high survival probability at 5 years [START_REF] Siegel | Cancer Statistics[END_REF]. However, 15% of the patients with early-stage breast cancer (eBC) will suffer from distant metastatic relapse after surgery, with limited treatment options [START_REF] Holleczek | Risk of Loco-Regional Recurrence and Distant Metastases of Patients with Invasive Breast Cancer up to Ten Years after Diagnosis -Results from a Registry-Based Study from Germany[END_REF][START_REF] Steeg | Targeting Metastasis[END_REF]. Prevention of metastatic relapse is the purpose of adjuvant (post-operative) systemic therapies designed to eradicate the minimal residual disease. Such therapies, which include chemotherapy, and/or hormone therapy in hormone receptor-positive tumor and/or trastuzumab in human epidermal growth factor 2 (HER2)-positive tumors [START_REF] Cardoso | Early Breast Cancer: ESMO Clinical Practice Guidelines for Diagnosis, Treatment and Follow-Up[END_REF], have substantially improved the metastasis-free and overall survivals [START_REF] Cameron | 11 Years' Follow-up of Trastuzumab after Adjuvant Chemotherapy in HER2-positive Early Breast Cancer: Final Analysis of the HERceptin Adjuvant (HERA) Trial[END_REF][START_REF] Ebctcg) | Comparisons between Different Polychemotherapy Regimens for Early Breast Cancer: Meta-Analyses of Long-Term Outcome among 100 000 Women in 123 Randomised Trials[END_REF][START_REF] Organisation | Controlled Trial of Tamoxifen as Adjuvant Agent in Management of Early Breast Cancer Interim Analysis at Four Years[END_REF]. Nevertheless, the clinical outcome of eBC patients is heterogeneous. Current routine prognostic features are mainly age, lymph node status, tumor size and grade, and HR and HER2 statuses.

However, several critical issues remain such as the identification of patients who would have been cured by surgery and radiotherapy alone, thus avoiding the use of toxic chemotherapy [START_REF] Pondé | Progress in Adjuvant Systemic Therapy for Breast Cancer[END_REF]. The current relapse risk assessment models are simple regressions based on the above-cited biological parameters (BP). Examples are the Nottingham Prognostic Index [START_REF] Galea | The Nottingham Prognostic Index in Primary Breast Cancer[END_REF] and the PREDICT score [START_REF] Wishart | A Population-Based Validation of the Prognostic Model PREDICT for Early Breast Cancer[END_REF][START_REF] Wishart | PREDICT: A New UK Prognostic Model That Predicts Survival Following Surgery for Invasive Breast Cancer[END_REF]. As of today, multiparameter genomic tests with elaborate gene expression signatures (e.g. MammaPrint [START_REF] Buyse | Validation and Clinical Utility of a 70-Gene Prognostic Signature for Women With Node-Negative Breast Cancer[END_REF][START_REF] Van 't Veer | Gene Expression Profiling Predicts Clinical Outcome of Breast Cancer[END_REF], Oncotype DX [START_REF] Sparano | Gene Expression Assays in Early-Stage Breast Cancer[END_REF] or Endopredict [START_REF] Gradishar | Predicting Expected Absolute Chemotherapy Treatment Benefit in Women with Early-Stage Breast Cancer Using a 12-Gene Expression Assay[END_REF]) can be used in clinical practice for predicting the clinical course. However, their use remains limited, in part due to their expensive price. More recently, machine learning algorithms have being developed for prognosis [START_REF] Kourou | Machine Learning Applications in Cancer Prognosis and Prediction[END_REF], but few of them focused on the prediction of breast cancer relapse [START_REF] Kim | Development of Novel Breast Cancer Recurrence Prediction Model Using Support Vector Machine[END_REF]. In addition, these approaches are agnostic and do not rely on biological knowledge.

In contrast, mechanistic models of metastatic development have been developed during the last decades, integrating the pathophysiology of the metastatic process [START_REF] Benzekry | Artificial Intelligence and Mechanistic Modeling for Clinical Decision Making in Oncology[END_REF]. They have been used to estimate the occult metastatic burden at diagnosis after the resection of the primary tumor [START_REF] Hanin | Suppression of Metastasis by Primary Tumor and Acceleration of Metastasis Following Primary Tumor Resection: A Natural Law?[END_REF][START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF][START_REF] Koscielny | A Simulation Model of the Natural History of Human Breast Cancer[END_REF], to predict the impact of individual treatments in pancreatic cancer [START_REF] Haeno | Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Optimum Treatment Strategies[END_REF] or to describe brain metastasis in lung cancer [START_REF] Bilous | Quantitative Mathematical Modeling of Clinical Brain Metastasis Dynamics in Non-Small Cell Lung Cancer[END_REF]. Following pre-clinical validation [START_REF] Baratchart | Computational Modelling of Metastasis Development in Renal Cell Carcinoma[END_REF][START_REF] Benzekry | Modeling Spontaneous Metastasis Following Surgery: An In Vivo-In Silico Approach[END_REF], we showed in a previous work [START_REF] Nicolò | Machine Learning and Mechanistic Modeling for Prediction of Metastatic Relapse in Early-Stage Breast Cancer[END_REF], that a mechanistic approach based on simulation of the metastatic disease could be used to create a predictive tool of breast cancer relapse. Building upon this work, we here tested our mechanistic models on three datasets of eBC patients who did not receive any adjuvant systemic therapy. This allowed us to calibrate the models using mixed-effects modeling from data of the natural history of the disease. We show -through a careful univariate analysis -that our model is able to describe the biological links between BP and processes of the metastatic disease. Then we propose a model selection procedure to establish the best covariate structures to use in prediction. Eventually, we establish the predictive performances of the selected models for each dataset.

Materials and methods

Patient datasets

The data consisted of distant metastasis-free survival (DMFS) information and main prognostic variables for patients with operated eBC from three databases. Inclusion criteria were: invasive breast carcinoma, early-stage, treated with primary surgery followed or not by adjuvant radiotherapy, without any adjuvant systemic therapy (hormone therapy, chemotherapy, trastuzumab), with clinicopathological data and follow-up available for DMFS. The patients who did not experience distant metastatic relapse were censored at the time of death or last follow-up.

The first dataset contained data from 591 women who were treated at the Bergonié institute (Bordeaux, France) between 1989 and 1993. The clinicopathological parameters were: age at the time of diagnosis, pathological tumor size, axillary lymph node status, tumor grade, and expression of estrogen (ER) and progesterone (PR) receptors, HER2 and Ki67, based on immunohistochemistry (IHC) assays. The tumors were considered ER-or PR-positive when more than 1% of the cells showed expression of the corresponding receptor in IHC, Ki67 high when 14% or more of the cells expressed the marker (and Ki67 low otherwise), and HER2-positive when the IHC score was 3+ or 2+ with 60% or more of the cells expressing HER2.

The second dataset included data from 676 patients extracted from our clinically annotated database (8, 982 invasive breast cancer samples) made from aggregation of 36 public gene expression datasets [START_REF] De Nonneville | NDRG1 in Aggressive Breast Cancer Progression and Brain Metastasis[END_REF]. This set included the same clinicopathological annotations as the Bergonié set.

The third dataset was composed of 167 patients treated between 1980 and 1990 at the public hospital of Marseille (AP-HM), France. Information on individual DMFS, age status, pathological tumor size, axillary lymph node status and grade were available with the same definition as for the other data sets. Protein dosage information, based on biochemical assays, was available for ER, PR, Urokinase Plasminogen Activator (UPA) and Plasminogen Activator Inhibitor-1 (PAI-1) as well as the enzymatic activity for Thymidine Kinase 1 (TK). Tumors were considered ER-or PR-positive when the quantity of respective proteins was greater than 15 femtomoles per milligram of protein.

Missing values could not be imputed uniformly across all three datasets and were thus removed from each dataset: 55 patients were removed from the Bergonié dataset (646 patients initially), 16 from the IPC dataset (692 patients initially), and 7 from the AP HM dataset (174 patients initially).

Mechanistic modeling of the metastatic process

The model has been introduced and extensively described previously [START_REF] Nicolò | Machine Learning and Mechanistic Modeling for Prediction of Metastatic Relapse in Early-Stage Breast Cancer[END_REF] but we include here a brief description for self-consistency. We consider a simple description of the natural history of eBC, starting at time 𝑡 = 0 with one cell. The primary tumor (PT) grows following a Gompertz law [START_REF] Norton | A Gompertzian Model of Human Breast Cancer Growth1[END_REF]:

𝑉 𝑝 (𝑡) = exp ( 𝛼 𝑏 (1 -𝑒 -𝑏𝑡 ) )
where 𝑉 𝑝 is the number of cells in the tumor at time t, 𝛼 is the specific growth rate (i.e. 𝑉 -1 𝑝 ⋅ 𝑑𝑉 𝑝 /𝑑𝑡, expressed in d -1 ) at 1 cell and 𝑏 is the exponential decay parameter of the initial growth rate (unitless) [START_REF] Vaghi | Population Modeling of Tumor Growth Curves and the Reduced Gompertz Model Improve Prediction of the Age of Experimental Tumors[END_REF]. These assumptions implicate that the tumor size converges to a theoretical limit 𝐾 = exp (𝛼𝑏) cells when 𝑡 → +∞. To avoid overparametrization and based on biological evidence [START_REF] Koscielny | A Simulation Model of the Natural History of Human Breast Cancer[END_REF][START_REF] Coumans | Detection of Cancer before Distant Metastasis[END_REF][START_REF] Klein | Parallel Progression of Primary Tumours and Metastases[END_REF][START_REF] Spratt | Decelerating Growth and Human Breast Cancer[END_REF], we fixed 𝐾 to 10 12 cells and considered 𝛼 as the only free parameter for growth (𝑏 being computed using the previous equation). At the time of diagnosis (𝑡 𝑑𝑖𝑎𝑔 ), the patient undergoes a surgery and the PT is removed.

During the course of the pre-surgical period, we assume that all cells from the primary tumor have an instantaneous probability of dissemination of 𝜇 (expressed in cell -1 d -1 ). The dissemination rate of the tumor is then:

𝑑(𝑉 𝑝 ) = 𝜇𝑉 𝑝 ,
leading to the following continuous expression for the total number of metastases at time 𝑡:

𝑁 𝑐𝑜𝑛𝑡 (𝑡) = ∫ 𝑡 0 𝑑 ( 𝑉 𝑝 (𝑠) ) 𝑑𝑠 = 𝜇 ∫ 𝑡 0 𝑉 𝑝 (𝑠) 𝑑𝑠
To be consistent with the biological reality, the numerical implementation considers that the number of metastases is an integer, given by:

𝑁(𝑡) = ⌊𝑁 𝑐𝑜𝑛𝑡 (𝑡)⌋,
where ⌊𝑥⌋ is the integer part of 𝑥. The metastases are assumed to also follow a Gompertz growth law with the same parameters as for the PT [START_REF] Benzekry | Modeling Spontaneous Metastasis Following Surgery: An In Vivo-In Silico Approach[END_REF].

To define the time to distant metastatic relapse (TTR), we considered that metastases are detected as soon as they reach a detectability threshold 𝑉 𝑑𝑒𝑡𝑒𝑐𝑡 taken to correspond to a tumor of 5 mm (detecting limit in imaging) [START_REF] Kundel | Predictive Value and Threshold Detectability of Lung Tumors[END_REF][START_REF] Macmahon | Guidelines for Management of Small Pulmonary Nodules Detected on CT Scans: A Statement from the Fleischner Society[END_REF]. From the size 𝑉 𝑑𝑖𝑎𝑔 of the PT at diagnosis and the growth parameters, we can compute 𝑡 𝑑𝑖𝑎𝑔 , the time between the initiation of the disease and the diagnosis and 𝜏 𝑣 𝑖𝑠, the time needed for a metastasis to reach 𝑉 𝑑𝑒𝑡𝑒𝑐𝑡 [START_REF] Nicolò | Machine Learning and Mechanistic Modeling for Prediction of Metastatic Relapse in Early-Stage Breast Cancer[END_REF].

Since the first metastasis emitted will be the first to reach the visibility threshold, the time to relapse is given as a function of 𝑉 𝑑𝑖𝑎𝑔 and the two mathematical parameters (MP) 𝛼 and 𝜇 (see Figure A where 𝜀 𝑖 ∼ 𝒩 (0, 𝜎 2 ) is the residual error with standard deviation 𝜎.

We also assumed a log-normal distribution of the individual MPs, in the population, with a linear effect of the covariates (BP at diagnosis), denoted by the vector 𝐶 𝑖 :

{ log 𝛼 𝑖 = log 𝛼 𝑝𝑜𝑝 + 𝛽 𝛼 ⋅ 𝐶 𝑖 + 𝜂 𝑖 𝛼 log 𝜇 𝑖 = log 𝜇 𝑝𝑜𝑝 + 𝛽 𝜇 ⋅ 𝐶 𝑖 + 𝜂 𝑖 𝜇
where 𝛼 𝑝𝑜𝑝 and 𝜇 𝑝𝑜𝑝 are the typical values of 𝛼 and 𝜇 in the population, 𝛽 𝛼 and 𝛽 𝜇 are the vectors of the covariate effects, and 𝜂 𝑖 = ( 𝜂 𝑖 𝛼 , 𝜂 𝑖 𝜇) are the random-effects, i.e., independent identically distributed random variables with distribution 𝒩 2 (0, Ω). The latter quantify inter-individual variability.

From the survival function implicitly defined by the structural error model and using a likelihood definition compatible with censored data, we defined a maximum likelihood estimator for the mixed-effect model (see previous work for technical details [START_REF] Nicolò | Machine Learning and Mechanistic Modeling for Prediction of Metastatic Relapse in Early-Stage Breast Cancer[END_REF]).

The relative standard-errors of the population-level MPs of the models -𝛼 𝑝𝑜𝑝 , 𝜇 𝑝 𝑜𝑝, Ω, and 𝜎 2 -were obtained by 100-replicates bootstrap and used to assess parametric identifiability.

Variable selection

The definition of the covariate effects allows for each BP to potentially influence the distribution of 𝛼, 𝜇, both, or none, depending on the corresponding coefficients 𝛽 𝛼 and 𝛽 𝜇 . To identify the impact of the BPs on the MPs, we used a two-step approach. First, we performed a univariate analysis in which we tested for significant effects in either 𝛼 or 𝜇, using models including only one covariate on one MP. Specifically, for a BP 𝐶 𝑘 , we tested for the null hypotheses 𝐻 0 ∶ 𝛽 𝑘,𝛼 = 0 or 𝛽 𝑘,𝜇 = 0. The univariate models were assessed in 100-samples bootstraps. The standard deviation of the bootstrap distributions was used to evaluate the precision of the MP estimation. For each covariate, we tested if the corresponding coefficient was significantly non null with a Wald test using the bootstrap estimate of the coefficient standard error.

In the second step, the optimal covariate model for each dataset was selected using a backward elimination procedure based on the Bayesian information criteria (BIC) adapted for the selection of covariates in mixed-effects models [START_REF] Delattre | A Note on BIC in Mixed-Effects Models[END_REF]. Specifically, we started from a model containing all statistically significant covariates for both 𝛼 and 𝜇. Then we iteratively generated all possible nested models with one covariate less and selected the model with the minimal BIC.

To verify that 𝛽 was correctly identifiable, the significance of each covariate coefficient in the selected multivariate model was then re-assessed in the multivariate model with a Wald test, based on 100-samples bootstrap estimation of the standard deviation.

Individual predictions

Individual predictions of survival curves Ŝ 𝑖 were obtained by taking the empirical expectation of the survival function with respect to the inter-individual variability over 𝑁 𝑠 𝑖𝑚

Prediction performance metrics

To assess the prediction performance of the models at a fixed time point, we used calibration plots. These were obtained by predicting the survival probabilities at landmark times 𝑡 𝑙 for all patients Ŝ 𝑖 (𝑡 𝑙 ), binning these into 8 quantiles groups and computing the median prediction in each bin. For each bin, this prediction was plotted against the actual data group DMFS at time 𝑡 𝑙 estimated by the Kaplan-Meier method [START_REF] Austin | Graphical Calibration Curves and the Integrated Calibration Index (ICI) for Survival Models[END_REF][START_REF] Kaplan | Nonparametric Estimation from Incomplete Observations[END_REF]. Performance of the models for prediction was also assessed by the concordance index using the individual predicted DMFS probability at 5 years to order the comparable pairs [START_REF] Harrell | Evaluating the Yield of Medical Tests[END_REF].

Numerical implementation of the model

The mechanistic model was implemented as an R package with high performance simulation code in C++. All model simulations used a time step of 20 days and were performed with tumor size expressed in number cells. Tumor diameter data were converted into numbers of cells assuming spherical shape and a cell density of 106 cells per mm3 [START_REF] Spratt | Decelerating Growth and Human Breast Cancer[END_REF][START_REF] Spratt | Rates of Growth of Human Solid Neoplasms: Part I[END_REF].

Parameter identification was performed using the stochastic approximation of expectation maximization algorithm [START_REF] Lavielle | Mixed Effects Models for the Population Approach[END_REF] implemented in the saemix R package version 3.0 [START_REF] Comets | Parameter Estimation in Nonlinear Mixed Effect Models Using Saemix, an R Implementation of the SAEM Algorithm[END_REF]. All computations were performed with R version 4.0.4 [START_REF] Team | A Language and Environment for Statistical Computing[END_REF].

Results

Biological parameters discrepancies between the datasets

The distributions of the BPs were different in the three datasets (Table 1). There was significantly more Ki67 high patients in the IPC dataset (chi-squared test, 𝑝 < 0.001) and more node-positive patients in the Bergonié dataset when compared to the IPC dataset (𝑝 < 0.001) and the AP-HM dataset (𝑝 < 0.001). The distribution of the grade values was significatively different between the Bergonié and IPC datasets (𝑝 < 0.001), with more low-grade tumors in the Bergonié data and more high-grade tumors in the IPC data. The proportion of ER-positive and PR-positive patients across the three datasets was also significatively different and no statistical difference was found for the distribution of HER2 in the Bergonié and IPC data.

The distribution of primary tumor size appeared log-normal for the three datasets, with a smaller median in the Bergonié dataset than in the IPC set (𝑝 < 0.001, Brown-Mood median test) whereas no statistical differences in median could be found between the IPC and AP-HM datasets (Supplementary Figure ??). The three datasets also exhibited differences in DMFS (Supplementary Figure ??), with lower DMFS in the IPC dataset than in the Bergonié dataset (Cox regression hazard ratio 𝐻𝑅 = 2.2, 𝑝 < 0.001) or in the AP-HM dataset (𝐻𝑅 = 1.8, 𝑝 < 0.001). No significant difference in DMFS was found between the Bergonié and AP-HM datasets (𝐻𝑅 = 1.3, 𝑝 = 0.157).

Inter-individual variability of the mathematical parameters accurately describes distant metastasis-free survival curves

We first used our mechanistic model without covariates to see if individual variability in the pathological tumor size at diagnosis (included as a direct parameter of the model), associated to log-normal inter-individual variability of the MPs 𝛼 and 𝜇, was able to describe the observed TTR in the three datasets [START_REF] Siegel | Cancer Statistics[END_REF]. Figures 1B,1C, 1D show the descriptive performances of the mechanistic models. For all three models, the model-based population survival curves correctly described the observed DMFS data and remained within the confidence intervals of the Kaplan-Meier estimators, except for a slight underestimation of the DMFS for small times.

The population values 𝑙𝑜𝑔 𝛼 𝑝𝑜𝑝 and 𝑙𝑜𝑔 𝜇 𝑝𝑜𝑝 were estimated with good accuracy. The relative standard error (RSE) was 10.7% on 𝑙𝑜𝑔 𝛼 𝑝𝑜𝑝 in the AP-HM dataset and the RSEs for the remaining parameters were below 10% in all three datasets (Table 2).

Mechanistic modelling yields biological insight on the impact of biological parameters on metastasis

To study the association of the BPs with either growth or dissemination, we tested all univariate models with effect of a BP as a covariate on either 𝛼 or 𝜇 (Table 2). On the AP-HM dataset, the TK concentration had a significant effect on 𝛼 and the concentrations of PAI 1 and UPA had a significant effect on 𝜇. On the Bergonié and IPC datasets, the ER, PR, HER2, and Ki67 statuses had a significant effect on 𝛼 and 𝜇, and the lymph node status had a significant effect on 𝜇. The grade of the tumor was also significantly associated with both 𝛼 and 𝜇 on the IPC dataset.

For each dataset, an optimal BP set was selected using a backward elimination procedure starting from the model including all significant BPs (Supplementary Table S1). Model selection was performed based on the BIC to compromise between model performances and number of covariates. For each dataset, the covariate model with the lowest BIC amongst the nested models was kept as the best model (Figure 2). Ties on the minimum (differences on the BIC < 4) were resolved by choosing the larger model, as the BIC is known to be favor smaller models when compared to other criteria [START_REF] Stoica | Model-Order Selection: A Review of Information Criterion Rules[END_REF]. The final models contained the ER and Ki67 statuses on 𝛼 for the Bergonié dataset (Figure 2A), a more complex model with effects of the grade and PR on 𝛼 and the HER2 and nodes statuses on 𝜇 for the IPC dataset (Figure 2B), and only PAI 1 on 𝜇 for the AP-HM dataset (Figure 2C).

Estimation of the population parameters for the selected models is presented in Table 3. The effects of PAI 1 on 𝜇 for the AP-HM data and grade 3 for the IPC dataset were estimated with good precision (𝑅𝑆𝐸 < 30%), whereas the effects of the other variables were estimated with larger but still acceptable uncertainty (𝑅𝑆𝐸 < 50%). Only the effect of the ER status on 𝛼 for the Bergonié data was estimated with high RSE (52%). In the final models, all BPs were significant in multivariate analysis as covariates (Wald test). The coefficients of all BPs were in the same order of magnitude as the corresponding inter-individual variability standard deviation, confirming that the estimated coefficients have significant impact on the individual parameter distributions.

Covariate models accurately describe DMFS curves

To assess the ability of the models to describe different subpopulations of patients in each dataset, we analyzed the model predicted DMFS for different subgroups of patients (Figure 3). On the Bergonié dataset, the model captured well the DMFS difference between the ER-positive and ER-negative patients and between the Ki67 high and Ki67 low patients (Figure 3A-3B). On the IPC dataset, the model was able to correctly describe the differences of DMFS when comparing patients with different grades, HER2, based on the presence of invaded lymph nodes, or between PR-positive and PR-negative patients (Figure 3C-F). For the AP-HM dataset, the model captured very well the DMFS for patient in the upper PAI 1 tercile, was adequate for the second tercile but seemed to slightly overestimate DMFS for larger times, and the DMFS for the lower tercile of PAI 1 was underestimated by the model (Figure 3G).

Predictive performances

Next, we evaluated the prediction performances of the models. Calibration curves indicated excellent individual predictive power for all three models at the 5-, 10-and 15years landmark times (Figure 4). Calibrations at 5 years (Figure 4A-C 4D-F) were good for all three models and for all DMFS probabilities, with a slight trend to overestimate DMFS probability at 10 years for the higher probability group. In the Bergonié dataset, the model slightly underestimated DMFS at 15 years for smaller probabilities while the over-estimation of higher probabilities was more pronounced (Figure 4G). Whereas for the IPC (4H) and AP-HM (Figure 4I) datasets, the calibration at 15 years was still good.

We then computed Harrell's concordance index (c-index) as another measure of the prediction performance [START_REF] Harrell | Evaluating the Yield of Medical Tests[END_REF]. The standard deviation for the c-index was computed in 10folds (8-folds for the AP-HM dataset) cross-validation. The performances were modest with a c-index of 0.63 (95% confidence interval (CI) [0.44-0.83]) in the Bergonié dataset, 0.71 (95% CI [0.42, 0.99]) in the AP-HM dataset and 0.60 (95% CI [0.54, 0.80]) in the IPC dataset.

Discussion

Classical statistical models of metastatic risk, although able to detect correlation between biomarkers and outcome, fail to give causal insights about the mechanisms at stake. Several genomics-based prognostic tools are commercially available (Oncotype DX Recurrence Score, Prosigna Risk of Recurrence score, EndoPredict, Breast Cancer Index) for estimation of the recurrence risk in HR-positive and HER2-negative eBC. However, the cost of these tests limits their clinical use. Our approach helps to provide mechanistic information from routine clinical markers. We used a simple mechanistic model of metastatic development based on two processes -growth and dissemination -to analyze DMFS data from three different datasets of eBC patients, two of which contained only routinely available data, whereas the third contained also non-routine markers (UPA, PAI-1 and TK), known to have biological roles in the metastatic process [START_REF] Jänicke | Randomized Adjuvant Chemotherapy Trial in High-Risk, Lymph Node-Negative Breast Cancer Patients Identified by Urokinase-Type Plasminogen Activator and Plasminogen Activator Inhibitor Type 1[END_REF]. We not only correctly described the DMFS in the three populations, but also showed that our method could be used to link biological features with specific parts of the metastatic process.

Studying the effect of the biological parameters on the population distribution of the growth parameter 𝛼 and the dissemination parameter 𝜇 allowed us to associate each predictive feature with one (or both) aspect(s) of the metastatic process. Specifically, based on the data from the AP-HM cohort, our model supports the association of protease UPA and its inhibitor PAI 1 with metastatic dissemination potential. This is consistent with previous pre-clinical studies [START_REF] Duffy | uPA and PAI-1 as Biomarkers in Breast Cancer: Validated for Clinical Use in Level-of-Evidence-1 Studies[END_REF]. The association of Thymidine kinase 1, an enzyme involved in DNA synthesis allowing cell division [START_REF] Sherley | Regulation of Human Thymidine Kinase during the Cell Cycle[END_REF] known to be associated with larger tumors at diagnosis in eBC [START_REF] Broët | Thymidine Kinase as a Proliferative Marker: Clinical Relevance in 1,692 Primary Breast Cancer Patients[END_REF], was found to impact the growth parameter 𝛼, again consistently with the biology. On the two other datasets (Bergonié and IPC), the model specifically associated the presence of invaded lymph nodes with 𝜇. Overall, these findings support the ability of our mechanistic model to identify the biological role of specific markers.

Based on the first step of our univariate analysis, we proposed a variable selection method to establish the best combination of biological parameters to include in a predictive model for each dataset. Using the BIC, we selected the best compromise between goodness-of-fit and the number of biological parameters. The selected models also improved the insights given by the univariate approach. In the Bergonié dataset, the model pinpointed the Ki67 marker as an important predictor of the growth parameter 𝛼 and as the least important predictor of the dissemination parameter 𝜇, in accordance with the established biological role of Ki67 as a marker of cell division and tumor proliferation [START_REF] Dowsett | Assessment of Ki67 in Breast Cancer: Recommendations from the International Ki67 in Breast Cancer Working Group[END_REF]. Similarly, in the IPC dataset, the lymph node status was selected in 𝜇 (and was the vari-able to be eliminated in 𝜇 in the Bergonié dataset). The hormone receptor statuses were relevant in 𝛼 for both the IPC (where PR status was selected) and the Bergonié (where ER status was selected) datasets. Their effect on 𝜇 was less clear, since no hormone receptor status was selected on 𝜇 but PR status persisted in 𝜇 up until late steps of the elimination process.

The prediction performances of the best model for each dataset were mitigated, with very good performances in calibration at various time points, but surprisingly low c-index values. The c-index computation may not be very accurate considering that the crossvalidation sampling was not stratified on the event indicator variable, meaning that the number of events and thus the number of comparable pairs for the c-index computation was probably different between the cross-validation folds.

The differences in DMFS within each dataset (Supplementary Figure 2), the measurement methods (IHC staining for Bergonié, mRNA expression for IPC, and protein dosage for AP-HM) or variable availability prevented the possibility to properly compare the results across datasets. We chose to study routine clinical and biological data to gather as many patients as possible and to match as closely as possible the information available in the clinic.

The prediction at the individual level still needs further investigation to be up-part with the existing agnostic models. In particular a larger and more homogenous dataset with more patients and raw values of the markers instead of dichotomized categories could give better results. The next step of the development of our model should be to integrate the impact of adjuvant treatment on the individual risk of metastatic recurrence. With a better identification of the individual values of the MP in a mechanistic model of the metastatic process under the course of adjuvant therapies would give a relevant framework to tackle the problem of identifying patients who could avoid cytotoxic chemotherapy, or limit its extent to a minimal number of cycles. For each dataset, the base models correspond to the case with no BP effect (aside from the pathological tumor size). The full models correspond to the best models from the selection procedure, where conditionally to the vector of covariates (𝐶 𝑖 ) included in the model, the mathematical parameters followed a log-normal distribution such that, log 𝛼 𝑖 and log 𝜇 𝑖 are gaussian with respective mean log 𝛼 𝑝𝑜𝑝 + 𝛽 𝛼 ⋅ 𝐶 𝑖 and log 𝜇 𝑝𝑜𝑝 + 𝛽 𝜇 ⋅ 𝐶 𝑖 (where 𝛽 𝛼 and 𝛽 𝜇 are vectors of the BP specific coefficients) and respective standard deviation 𝜔 𝛼 and 𝜔 𝜇 . For all models, the log-residual error on time to distant metastatic relapse was assumed to follow a centered gaussian distribution with variance 𝜎 2 . Estimation was performed using the stochastic approximation of expectation maximization algorithm. Relative standard errors (RSE) were computed from a 100 replicates bootstrap.

Bergonié

BIGARRÉFigure 1 .Figure 2 .

 12 Figure 1. Mechanistic model of the time metastatic relapse.A. Overview of the mechanistic model. The model prediction of the time to metastatic relapse (TTR) is computed from the size of the primary tumor at diagnosis and two mathematical mechanistic parameters, α controlling the growth rate of the primary tumor and metastases, and μ controlling the seeding of new metastases. The scheme shows a unitless simulation of the model. B -D. Fits of the mechanistic model (without covariate effects) on the Bergonié (B), IPC (C) and AP-HM (D) datasets. Each panel presents for one dataset, the DMFS and the model's prediction of the survival function in the population. The model was trained and evaluated on the full dataset.

Figure 3 .Figure 4 .

 34 Figure 3. Model predictions in stratified groups. Group-comparison of the distant metastasis-free survival data (Kaplan-Meier estimate, solid line and 95% confidence interval, colored band) and the model mean prediction of the metastasis-free survival function. Bergonié dataset, patients stratified on the ER (A) or Ki67 (B) status. Paoli-Calmettes Institutes (IPC) dataset, patients stratified by grade (C), HER2 (D) status, invaded lymph nodes (E) or PR (F) status. G. APHM dataset, patients stratified by PAI-1 tertiles.

Statistical mixed-effects population model Given

  ). individual values 𝛼 𝑖 and 𝜇 𝑖 for the 𝑖-th patient with observed size of the PT at diag-

	relapse 𝑇 𝑖 , to ensure positivity,
	log(𝑇 𝑖 ) = log ( 𝑇 𝑇 𝑅(𝛼 𝑖 , 𝜇 𝑖 ; 𝑉 𝑖 𝑑𝑖𝑎𝑔 ) ) + 𝜀 𝑖 ,
	𝑇 𝑇 𝑅(𝛼, 𝜇, 𝑉 𝑑𝑖𝑎𝑔 𝑠) =	⎧ ⎪ 𝜏 𝑣𝑖𝑠 + arg min 𝑡 {𝑁(𝑡) ≥ 1} -𝑡 𝑑𝑖𝑎𝑔 , when 𝑁(𝑡 𝑑𝑖𝑎𝑔 ) ≥ 1 ⎨ ⎪ ⎩ +∞, otherwise
	nosis 𝑣 𝑖 𝑑𝑖𝑎𝑔 , we assumed a log-normal observation error model for the time to metastatic

  with 𝛼 𝑖 𝑗 and 𝜇 𝑖 𝑗 sampled from the distribution of 𝛼 𝑖 and 𝜇 𝑖 .

	Ŝ 𝑖 (𝑡) =	1 𝑁 𝑠𝑖𝑚	𝑁 𝑠𝑖𝑚 ∑ 𝑗=1	𝑆 ( 𝑡| 𝛼 𝑖 𝑗 , 𝜇 𝑖 𝑗 ; 𝑉 𝑖 𝑑𝑖𝑎𝑔 ) ,

BIGARRÉ et al. 2023 | Mechanistic modeling of metastatic relapse in early breast cancer | 5 of 21

replicates,

Table 1 .

 1 , N = 591 IPC, N = 676 APHM, N = 167 p-value Patient and disease characteristics. For categorical variables (age, presence of invaded lymph nodes, ER, PR , HER2 and Ki67 statuses, grade), number of patients in each dataset (and proportion of the dataset's values (%)). For continuous variables (pathological tumor size and PAI-1), median value (first quartile -third quartile) in each dataset. P-values correspond to the adequate test of identical distributions in all datasets with available information (chi-squared test for categorical variables, Kruskal-Wallis rank sum test for continuous variables).

	Pathological tumor size	15 (12, 20)	20 (15, 25)	20 (15, 26)
	Age			<0.001
	<50	147 (25%)	338 (71%)	35 (21%)
	≥50	444 (75%)	138 (29%)	132 (79%)
	Unknown	0	200	0
	Pathological lymph nodes			<0.001
	nodes-	364 (62%)	648 (96%)	147 (88%)
	nodes+	227 (38%)	28 (4.1%)	20 (12%)
	ER status			0.038
	ER-	134 (23%)	140 (21%)	50 (30%)
	ER+	457 (77%)	536 (79%)	117 (70%)
	PR status			<0.001
	PR-	189 (32%)	271 (40%)	80 (48%)
	PR+	402 (68%)	405 (60%)	87 (52%)
	HER2 status			0.7
	HER2-	523 (88%)	594 (88%)
	HER2+	68 (12%)	82 (12%)
	Grade			<0.001
	1	182 (31%)	131 (19%)
	2	267 (45%)	298 (44%)
	3	142 (24%)	247 (37%)
	Ki67 status			<0.001
	low	390 (66%)	273 (40%)
	high	201 (34%)	403 (60%)
	TK			60 (31, 154)
	UPA			0.79 (0.38, 1.33)
	PAI-1			3.4 (2.0, 5.3)

Table 2 . Values of the Parameters.

 2 The mathematical mechanistic parameters (𝛼 and 𝜇) were assumed to follow a log-normal distribution such that log 𝛼 𝑖 and log 𝜇 𝑖 are gaussian with respective mean log 𝛼 𝑝𝑜𝑝 and log 𝜇 𝑝𝑜𝑝 , and respective standard deviation 𝜔 𝛼 and 𝜔 𝜇 .

Data availability

The data from Institut Bergonié and AP-HM analyzed in this study are not publicly available due to patient privacy requirements but are available upon reasonable request from the corresponding author. The data from IPC is available online and has been previously described [START_REF] De Nonneville | NDRG1 in Aggressive Breast Cancer Progression and Brain Metastasis[END_REF].

Appendix 1 Supplementary Tables

The univariate models assumed a log-normal distribution independent of the covariate for one of the two computational biomarkers (𝛼 or 𝜇), and a log-normal conditional distribution with respect for the covariate for the other computational biomarker, with a median equal to the sum of a typical population value, and of the weighted covariate value. The estimation of the weight coefficient (beta), the relative standard error (R.S.E.) obtained by bootstrap with 100 repetitions, as well as the p-values for the corresponding Wald test, are presented in the table for all possible univariate models.

Bergonié


Yacine Bechane

Benoît Fiorina

Numerical analysis of turbulent flame enhancement by nanosecond repetitively pulsed plasma discharges

Keywords: Plasma-assisted combustion, Nanosecond Repetitively Pulsed discharges, Turbulent premixed flame, Lean flame enhancement, Large Eddy Simulation

Experimental studies of plasma-assisted combustion have shown that Nanosecond Repetitively Pulsed (NRP) discharges are a very efficient way to improve flame stability in lean regimes. These discharges generate a nonequilibrium plasma that induces a local heating and an important production of active species sufficient to enhance the combustion. The aim of this article is to understand the physical and chemical mechanisms involved in plasma-assisted turbulent combustion. High-performance computations of a lean bluff-body turbulent premixed methane-air flame, experimented at EM2C laboratory, are conducted for that purpose. Simulations are performed by combining a semi-empirical plasma model with an analytical combustion mechanism. Without electrical discharges, only a weak reaction zone behind the bluff body is observed. However, when the flame is assisted with the plasma, the flame power and surface increase significantly. A chemical analysis performed at the flame basis shows that the heat and the radical O produced by NRP discharges induce the dissociation of main burnt gases species H2O and CO2 into radicals OH, H, CO, O and species H2, O2. The species with long lifespan OH, H2 and O2 are convected from the center of the bluff-body recirculation zone to the flame front where they are consumed increasing the chemical reactivity. This causes a local increase of heat release rate that attaches the turbulent flame front which can develop downstream.

1 1. Introduction
The low temperature of lean flames prevents the NOx production but may induce incomplete combustion, instabilities and extinctions [START_REF] Candel | Combustion dynamics and control: Progress and challenges[END_REF].

An emerging solution to enable flame stabilization, suitable to a wide range of combustion applications, is to locally generate a plasma [START_REF] Criner | Stabilisation of turbulent lifted jet flame assisted by pulsed high voltage discharge[END_REF][START_REF] Kim | Flame stabilisation enhancement and nox production using ultra short repetitively pulsed plasma discharges[END_REF]. Among the various possible types of discharges to enhance combustion [START_REF] Starikovskaia | Plasma assisted ignition and combustion[END_REF][START_REF] Ju | Plasma assisted combustion: dynamics and chemistry[END_REF], the Nanosecond Repetitively Pulsed (NRP) discharges [START_REF] Pilla | Etude experimentale de la stabilisation de flammes propane-air de pre-melange par decharges nanosecondes impulsionnelles repetitives[END_REF] succeeded to stabilize ultra-lean premixed flames on several laboratory-scale hydrocarbon-air flames [START_REF] Pilla | Etude experimentale de la stabilisation de flammes propane-air de pre-melange par decharges nanosecondes impulsionnelles repetitives[END_REF][START_REF] Bak | Studies on the stability limit extension of premixed and jet diffusion flames of methane, ethane, and propane using nanosecond repetitive pulsed discharge plasmas[END_REF] but also in test rigs representative of a gas turbine environment at atmospheric [START_REF] Barbosa | Influence of nanosecond repetitively pulsed discharges on the stability of a swirled propane/air burner representative of an aeronautical combustor[END_REF][START_REF] Vignat | Improvement of lean blow out performance of spray and premixed swirled flames using nanosecond repetitively pulsed discharges[END_REF] and high pressure [START_REF] Sabatino | Enhancement of the lean stability and blow-off limits of methane-air swirl flames at elevated pressures by nanosecond repetitively pulsed discharges[END_REF][START_REF] Sabatino | Actuation efficiency of nanosecond repetitively pulsed discharges for plasmaassisted swirl flames at pressures up to 3 bar[END_REF]. This technique is very efficient because the energy consumption of the plasma is typically less than 1% of the power released by the flame.

The principle of flame stabilization by NRP discharges is illustrated in Fig. 1. NRP discharges are generated by applying a high voltage pulse (about 10kV) during approximatly 10 ns between two electrodes at a repetition frequency of the order of 10kHz. The NRP discharges cause first, at the nanosecond scale, an ultrafast increase of gas temperature and species dissociation [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF]. This is followed by a slower gas heating process, appearing at the millisecond scale, induced by the relaxation of the nitrogen molecule's vibrational states. It has been shown experimentally [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF] and numerically [START_REF] Popov | Fast gas heating initiated by pulsed nanosecond discharge in atmospheric pressure air[END_REF] that about 35% of the discharge energy dissociates up to 50% of the oxygen molecules contained in the inter-electrode region in O atoms. While it has been observed empirically that the NRP discharge is efficient when located in small gases recirculation zone, the fundamental mechanisms of plasma-assisted flame stabilization are not yet well understood. In particular, the physical and chemical processes that relate the O atom production to the heat release enhancement in a turbulent flame has not been completely identified.

High-performance computing is an attractive way to understand and explain such complex fundamental mechanisms. Very few CFD simulations have been however carried out to understand the influence of NRP discharges on the flame structure. Bak et al. [START_REF] Bak | Plasmaassisted stabilization of laminar premixed methane/air flames around the lean flammability limit[END_REF] have simulated the plasma-assisted stabilization of a laminar methane/air flame by simultaneously solving the equations describing the discharge, the plasma, the flow, and the combustion kinetics. A similar computational strategy has been followed by Tholin et al. [START_REF] Tholin | Influence of fast-heating processes and o atom production by a nanosecond spark discharge on the ignition of a lean h2-air premixed flame[END_REF] to simulate the ignition of preheated quiescent H2/air mixtures excited by NRP discharges. While giving interesting insights into the impact of plasma on combustion, detailed simulations have been limited to 1-D and 2-D simulations in laminar flow conditions, because of their high computational cost. Computational domains have been also restricted to the plasma discharge vicinity, which is not sufficient to capture the large scales involved in a turbulent combustion chamber. Detailed simulations of plasma kinetics is today still too expensive to be directly in-cluded in 3-D turbulent flame simulations.

An alternative has been proposed by Castela et al. [START_REF] Castela | Modelling the impact of nonequilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows[END_REF] to reduce the computational cost of the description of the complex plasma kinetic processes by using a phenomenological model which captures the main impacts of NRP on the combustion: the local increase of temperature and the oxygen dissociation [START_REF] Castela | Modelling the impact of nonequilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows[END_REF]. Only one balance equation for the vibrational energy is added to the standard reactive flow balance equations, therefore the increase of CPU cost is limited. This approach has been succesfully applied to compute the ignition of flame kernel by NRP discharges [START_REF] Castela | A 3-d dns and experimental study of the effect of the recirculating flow pattern inside a reactive kernel produced by nanosecond plasma discharges in a methane-air mixture[END_REF]. This model has also been recently implemented in a Large Eddy Simulation (LES) flow solver to study the ignition sequence of a flame by a series of NRP discharges [START_REF] Bechane | Numerical investigations of turbulent premixed flame ignition by a series of nanosecond repetitively pulsed discharges[END_REF].

The objectives of this study are to give new insights into the plasma flame interaction mechanisms through the numerical analysis of the LES of plasma-assisted turbulent combustion. Section 2 presents the phenomenological plasma-assisted turbulent combustion model adapted with the LES formalism, while Section 3 gives details on the experimental configuration, called Mini-PAC, retained for performing the numerical investigation. Results are presented and analyzed in Sec.4 .

Modeling plasma-assisted turbulent combustion

Nanoseconds pulsed discharges model

The experimental [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF] and numerical [START_REF] Popov | Fast gas heating initiated by pulsed nanosecond discharge in atmospheric pressure air[END_REF] studies of NRP discharges applied in pure air or lean fresh hydrocarbon-air mixtures showed that the discharge pulses mainly induce electronic and vibrational excitation of nitrogen molecules N2. The nanosecond scale relaxation of the produced electronic states N2(A,B,C,a', . . . ) by dissociative quenching with O2 molecules leads leads to an ultra-fast increase of gas temperature and O2 dissociation into atomic oxygen [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF].The Relaxation of vibrational states of N2 occurs on a much longer times scales, milliseconds scales, and results in slow gas heating. The recent measurement performed by Minesi et al. [START_REF] Minesi | Plasma-assisted combustion with nanosecond discharges. i: Discharge effects characterization in the burnt gases of a lean flame plasma[END_REF] showed that when applied in a mixture of lean flame burnt gases the NRP discharges induce similar effects.

Based on these observations, Castela et al. [START_REF] Castela | Modelling the impact of nonequilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows[END_REF] developed a semi-empirical model that captures the effects of NRP dischrages on combustion without solving the electric field and the detailed plasma chemistry. Ultra-fast and slow processes are modeled by splitting into three contributions the rate of energy Ėp deposited at each electrical discharge :

Ėp = Ėp chem + Ėp heat + Ėp vib (1)
The chemical Ėp chem and thermal Ėp heat contributions refer to the amount of energy transferred into the ultra-fast dissociation of species and gas heating, respectively. Ėp vib is the rate of energy transferred to vibrational state N2. The model considers that only O2 is dissociated by the discharge into atomic O. The energy contributions are modeled as [START_REF] Castela | Modelling the impact of nonequilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows[END_REF]:

Ėp chem = η YO 2 Y f O 2 (1 - eO 2 eO ) Ėp (2)
Ėp heat = [α -η YO 2 Y f O 2 (1 - eO 2 eO )] Ėp (3)
Ėvib = (1 -α) Ėp (4)
where Y f O 2 is the mass fraction of O2 in the fresh mixture. e k is the internal energy of species in mass units. α is the fraction of pulse energy that induces the electronic excitation of N2 and subsequently the ultra-fast gas heating and ultra fast O2 dissociation. η is the fraction of pulse energy that leads to O2 dissociation, αη the fraction that goes into ultra-fast heating and finally (1 -α) the remainings fraction that goes into the slow heating.

According to [START_REF] Minesi | Plasma-assisted combustion with nanosecond discharges. i: Discharge effects characterization in the burnt gases of a lean flame plasma[END_REF], the fraction of energy α responsible of the ultra-fast heating and the ultra-fast dissociation of species when NRP discharges are applied in lean burnt gases is similar to what was obtained for NRP discharges in pure air [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF]. The fraction of energy that goes into ultra-fast heating (α -η) is also the same and is equal to 20%. Therefore, the parameters are here set to α = 55% and η = 35%.

LES balance equations

The LES reactive flow governing equations are combined with the Thickened Flame model for LES (TFLES) [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF], which artificially thicken the reactive layer to avoid under-resolution issues of the flame front. In addition to the mass and momentum equations, the following balance equations are solved for total energy e, vibrational energy e vib and species mass fractions Y k :

∂ρ e ∂t + ∂(ρ ui e) ∂xi = ∂ ∂xi Cp F Ξ µ Pr + (1 -S) µt P t r ∂ T ∂xi + ∂ ∂xi N k=1 F Ξ µ S c,k + (1 -S) µt S t c,k ∂ Y k ∂xi h s,k + ∂(σijui) ∂xi + Ėp chem + Ėp heat + Ṙp V T (5
)

∂ρ e vib ∂t + ∂(ρ ui e vib ) ∂xi = ∂ ∂xi µ S c,evib + µt S t c,evib ∂ e vib ∂xi + Ėp vib - Ṙp V T (6) ∂ρ Y k ∂t + ∂(ρ ui Y k ) ∂xi = ∂ ∂xi F Ξ µ S c,k + (1 -S) µt S t c,k ∂ Y k ∂xi + Ξ F W k ωc k + W k ωp k (7)
where φ and φ are Reynolds and Favre filtering of variable φ. µ, Sc and Pr are the viscosity, the Schmidt and Prandtl numbers, respectively, where t superscript denotes for turbulent quantities. The dynamical formulation of the thickening factor reads F = 1 + (Fmax -1)S, where the flame sensor S equals 0 and 1 outside an inside the reactive layer, respectively. Ξ is the sub-grid scale flame wrinkling. ωc k , the chemical reaction rate due to combustion, is modeled by a detailed mechanism. ωp k , the filtered chemical reaction rate due to plasma chemistry is modeled by assuming that the dissociation of O2 into O dominates:

ωp O = η YO 2 Y f O 2 Ėp eO WO (8)
ωp O2 = - WO WO2 ωp O ( 9
)
ωp k = 0 if k ̸ = O2, O (10)
As the plasma chemistry are much faster than flow time scales, sub-grid scale interactions between turbulence and plasma have been neglected. The term Ṙp V T in both Eq. ( 5) and Eq. ( 6) refers to the relaxation rate of the vibrational energy into gas heating and is closed using a Landau-Teller harmonic oscillator as indicated in [START_REF] Bechane | Numerical investigations of turbulent premixed flame ignition by a series of nanosecond repetitively pulsed discharges[END_REF].

The filtered rate of energy, Ėp chem , Ėp heat and Ėp vib are closed from Eq. ( 2), Eq. (3) and Eq. ( 4) by neglecting sub-grid scale contribution, i.e. by assuming that Ėp y = Ėp y ( YO 2 , Ėp ).

MiniPAC configuration

Experimental set-up

The MiniPAC configuration is a premixed methane-air bluff-body stabilized turbulent flame. The 200-mm-length injection tube has an inner diameter of 16 mm and the bluff-body diameter is 10 mm. Nanosecond discharges are produced downstream of the bluff-body between a horizontal tungsten anode of 1 mm diameter centered 5 mm above the grounded cathode which consists of a 2 mm height, 1 mm diameter pin inserted in the middle of the bluff-body.

The burner is fed with a methane-air mixture (ϕ = 0.8), injected at a constant volumetric flow rate ṁ = 17.4 m 3 h -1 which corresponds to a bulk velocity of U b = 43.3 m.s -1 and a Reynolds number of Re = 1.6 10 4 .

The discharges properties have been experimentally characterized [START_REF] Blanchard | Experimental and numerical characterization of a lean turbulent premixed flame stabilized by nanosecond discharges[END_REF]. They have a quasi-cylindrical shape of 5 mm length corresponding to the interelectrode distance, and a radius equal to 1.0 mm, a pulse energy of 2 mJ deposited in 50 ns and repeated at a frequency of 20 kHz. Photographies of non-assisted and plasma-assisted flames are shown in Fig. 2. Without discharge (left), a weak and incomplete combustion is restricted to the bluff-body wake. However, the application of NRP discharges allows the flame enhancement and stabilization (right).

Numerical set up

LES of the MiniPAC stabilization by discharges are performed using the YALES2 unstructured finite-volume low-Mach number solver [START_REF] Moureau | Design of a massively parallel cfd code for complex geometries[END_REF]. The time integration is performed using a fourth-order temporal scheme with a centred fourth order scheme for spatial discretization. The sub-grid Reynolds stresses tensor is closed using the Dynamic Smagorinsky model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. The analytical LU19 kinetic scheme involving 19 species and 184 reactions is used for the methane-air mixture chemistry [START_REF] Lu | A criterion based on computational singular perturbation for the identification of quasi steady state species: a reduced mechanism for methane oxidation with no chemistry[END_REF], species equations are coupled with turbulence using the Thickened Flame model for LES (TFLES) [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF] and the Charlette sub-grid scale flame wrinkling model [START_REF] Charlette | A power-law wrinkling model for les of premixed turbulent combustion: Part i non dynamic formulation and initial tests[END_REF]. The mesh consists of 140 millions tetrahedral elements refined in the plasma discharge zone, with a cells size equals 0.1 mm, to ensure that the small turbulence scales and laminar flame thickness are well resolved and no sub-grid scale model or TFLES coupling is necessary [START_REF] Bechane | Numerical investigations of turbulent premixed flame ignition by a series of nanosecond repetitively pulsed discharges[END_REF]. In the flame area, the cells size varies from 0.1 mm to 0.2 mm. The TFLES model is used in this region to capture the turbulence-flame interaction, triggered by a dynamic sensor.

The discharge properties have been set to be consistent with experimental measurements. It has a cylindrical shape filtered with the spatial function F (r) = erf c[(r/a) b ] where r is the radial distance from the discharge axis. a and b are geometric parameters, fitted to well reproduce the discharge radius and have a sufficient number of point at discharge interface to avoid high gradients that can cause numerical issues.

A reactive LES is first conducted without NRP discharges until reaching a steady-state solution. Then, from the same initial condition,which define t=0, two LES are performed: one still without activating plasma (case NP) and a second one with applying an uninterrupted series of NRP discharges in the burnt gases of the flame above the bluff-body at a frequency of 20kHz (case P).

Figure 3 shows the 2D mean axial velocity on the centerline plane obtained from case NP. An inner recirculation zone is clearly identified behind the bluffbody. The same analysis conducted on case P shows that the shape and the size of the recirculation are not significantly affected by the electric discharges.

Analysis of the plasma-assisted combustion mechanism

Quantitative comparisons between experiments and simulations have been performed to validate the modeling approach consistency. The comparisons of the gas temperature evolution at the center of the discharge in [START_REF] Blanchard | Experimental and numerical characterization of a lean turbulent premixed flame stabilized by nanosecond discharges[END_REF] and OH profile along the line passing through the center of the discharge (4.5 mm above the bluff body) [START_REF] Minesi | Thermal spark formation and plasmaassisted combustion by nanosecond repetitive discharges[END_REF] during the pulsing sequence show a very good agreement between the experiment and the numerical simulation with less than 10% difference when steady state is reached. The local effects of the NRP discharges applied on lean burnt gases are well capture by the model.

Global scale comparisons have also been performed.The direct experimental visualization of the OH* chemiluminescence against the simulated OH mass fraction conducted by Blanchard et al. show that the influence of plasma on the flame shape and length is also well predicted by the LES.

To understand the local and global interaction mechanisms between the plasma and the flame involved in plasma-assisted flame enhancement process a detailed numerical analysis is performed.

Global impact of NRP discharges

The global impact of NRP discharges on the total heat release is evidenced in Fig. 4, which plots the temporal evolution of the thermal flame power P = V (Ξ/F ) ωH dV , where ωH is the filtered heat release rate and V denotes the volumetric integral over the whole computational domain V . Without NRP discharges, the flame power oscillates around 700 W. However, when the electrical discharges are turned on, the total heat released by the combustion increases gradually during a transient period of approximately 10 ms (200 pulses) and then oscillates around a value of 1600 W.

To identify the reasons of the flame power increase, the LES-resolved R and total A flame surfaces are computed as [START_REF] Veynante | Turbulent combustion modeling[END_REF]:

R = V |∇c|dV and A = V Ξ|∇c|dV (11)
where c is a progress variable defined from the temperature to reach 0 and 1 in the fresh and burnt gases, respectively. Two observations are made from the analysis of the temporal evolution of R and A plotted in Fig. 5. First, R is very close to A and in practice 95% of the flame wrinkling is resolved at the LES scale. The grid resolution is here sufficiently fine to limit turbulent combustion modeling errors. Second, the flame surface evolution follows the same trend than the flame power plotted in Fig. 4.

The mean gains in power (ηP ) and total flame surface (ηA) induced by the NRP discharges are defined as: ηP = ⟨PP ⟩ / ⟨PNP ⟩ and ηA = ⟨AP ⟩ / ⟨ANP ⟩ (12) where ⟨φ⟩ denotes temporal averaging of variable φ. P and N P subscripts refer to simulations with and without plasma, respectively. To remove the transient state following by the discharges activation, the mean gains are computed for the period [7ms,20ms]. It comes that ηP = 2.13 and ηA = 1.54. The ratio ηA/ηP = 0.72, which means that approximately 70 % of the gain in flame power is explained by an increase of the total flame surface. As it will be shown further, the rest of the gain (30%) is explained by a local increase of the heat release (and the flame consumption speed) in the periphery of the bluff-body recirculation zone.

Impact of NRP discharges on the heat release field

Figure 6 shows two instantaneous snapshots from cases NP (left) and P (right), taken at the same time instant t= 20 ms, of the flame surface ( defined as the iso-surface of c * = 0.5) colored by the heat release. Without NRP discharges, only a weak reaction zone limited to the recirculation zone behind the bluff body is visible. Indeed, because of the low burning velocity under lean conditions, the flame front cannot propagate outside the recirculation zone. Figure 6 (right) plots the flame surface after 400 pulses. In comparison with case NP, the heat released by the flame increased significantly and the reaction zone extended downstream. As expected from the previous global analysis, the flame surface is developed. The pulsation of the electrical discharges enables the recirculation zone to act as a flame holder. A chemical analysis is conducted to understand how the NRP discharges, located at the center of the bluff-body, will attach the flame around the recirculation zone.

Chemical influence of NRP discharges in the recirculation zone

Average thermo-chemical quantities are computed over the period [10ms,20ms], for both cases NP and P. Radial profiles then are extracted at the position indicated in Fig. 3, located 4.5 mm above the bluff body, which crosses both the center of the discharges and the bluff-body recirculation zone. Mean radial heat release and methane consumption rates are plotted in Figs. 7 (top) and (bottom), respectively. The peak of heat release observed in both NP and P cases corresponds to the flame reactive layer, whose position was also evidenced in Fig. 6. The application of NRP discharges at the center of the bluff-body (x=0) causes a significant increase of heat release rate in the reactive layer located at 4 mm < r < 5 mm, which is accompanied by an increase of the methane consumption. The methane source term and the heat release rate are affected by the same factor (≈2.2), which means that the increase of the heat release rate is directly related to the increase of the methane consumption when NRP discharges are applied. Figure 8 shows the mean radial profile of temperature across the recirculation zone. For case NP, the temperature of the burnt gases contained within the recirculation zone is around 1950K which corresponds to the adiabatic flame temperature of a methane-air mixture at an equivalence ratio of 0.8. However, for the plasma-assisted flame, the mean temperature increases up to 3500K in the plasma zone, where the thermal energy is deposited. This energy is convected by the flow from the center to the periphery of the recirculation zone and induces a temperature increase in the flame front, i.e. at 4 mm < r < 5 mm.

Figures 9 and 10 plot the mass fraction of species H2O-CO2-O2-CO and OH-O-H2-H, repectively. Without plasma, the chemical composition within the recirculation zone corresponds to thermochemical equilibrium condition of the methane-air mixture with equivalence ratio of 0.8. When plasma is activated, both H2O and CO2 are highly consumed in the discharge zone. At the opposite, the mass fractions of the other species, in particular O and OH, increases significantly.

Figure 11 (left) and (right) shows the chemical pathway diagram in the plasma region during and after a pulse, respectively. The thickness of the arrow is proportional to the reaction rates integrated within the discharge zone. As discussed in [START_REF] Rusterholtz | Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges[END_REF][START_REF] Castela | Modelling the impact of nonequilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows[END_REF], O2 is mainly dissociated into O atoms during an electric pulse (Fig. 11 (left)), by the relaxation of nitrogen electronic state. The, the atomic oxygen produced by the discharge reacts with H2O contained in the burnt gases to produce OH and H. As shown in Figure 11 (right), chemical pathways are completely different after the pulse. CO2, the main reactant, is consumed by heat and radicals produced during the pulse, leading to a production of CO, H2, O2, O and OH. Fig. 9 and 10 show that, among the species produced in the discharge zone, those with a long lifespan such as OH, O2 and H2 are convected up to the flame front in large quantities. At the opposite, species with a shorter lifespan such as O, H does not survive outside the discharge zone.

A second chemical analysis, shown in Fig. 12, is performed in the flame reactive layer ( 4 mm < r < 5 mm) front to identify the CH4 consumption pathways without and with NRP discharges. Without plasma, CH4 is mainly consumed by OH, H and O trough the following three elementary reactions:

CH4 + OH → CH3 + H2O (R1) CH4 + H → CH3 + H2 (R2) CH4 + O → CH3 + OH (R3)
The impact of the plasma on the rate of progress of reactions R1, R2 and R3 is indicated in Fig. 12. As discussed previously, NRP discharges supply the reaction zone with OH, H2 and O2. The increase of OH concentration will enhance directly the reaction R1 by 51%. H2 and O2 are dissociated into H and O radicals in the flame reaction zone, increasing therefore also the rate of reactions R2 and R3 by 27% and 21%, respectively.

In conclusion, the discharges produces heat and species (OH,H2,O2), which are convected within the recirculation zone. Because of its longer lifespan, OH is the main plasma-induced species that reaches the flame front and enhance significantly the chemical reactivity at the flame basis surrounding the recirculation zone. Such local increase of the reactivity enhances the heat release in the periphery of the recirculation zone but also acts as a flame-holder which enables the turbulent flame surface development, increasing the total combustion efficiency. The results of the analysis are consistent with the experimental images in [START_REF] Xu | Thermal and hydrodynamic effects of nanosecond discharges in air and application to plasmaassisted combustion[END_REF] that show an enhancement of the flame when the OH produced by the discharge reaches the flame front.

Conclusion

The numerical analysis of Large Eddy Simulation of the Mini-PAC burner concludes to the following mechanisms of turbulent premixed flame stabilization by nanosecond repetitively pulsed (NRP) discharges:

1. The electrical pulses generate a plasma in the discharge zone, which produces O atom and heat. 2. During the pulse, O reacts with H2O to form mainly OH, through the reaction O+ H2O → 2OH. After the pulse, CO2 reacts in the interelectrode region to form H2 and O2. 3. Because of their long lifespan, OH, H2 and O2 are the main plasma-induced species convected to the flame front thanks to the bluff-body recirculation zone. 4. OH, H2 and O2 are consumed with CH4 in the flame reactive layer surrounding the recirculation zone and locally increase the flame heat release rate. This local increase of the heat release in the vicinity of the recirculation zone is responsible of 30% of the gain in flame power induced by the plasma. 5. The local increase of heat release rate anchored the turbulent flame front which can develop downstream. 6. The total flame surface is increased, increasing the combustion efficiency by 70%. This analysis is based on the assumption that O2 dissociation is the key chemical effect of plasma on combustion chemistry. Future works will aim to challenge this hypothesis by extending the plasma combustion model to account for the dissociation of other species present in the recirculating burnt gases, such as CO2 and H2O dissociation induced by the discharge plasma chemistry.

72 Chapter 3 -Figure 3 . 3 :Figure 3 . 3 :Fig. 1 :

 333331 Figure 3.3: Half portion of the TSF-A computational domain. A longitudinal cut of the mesh is shown and colored by the filtered progress variable reaction rate e !Yc .

Fig. 2 :

 2 Fig. 2: Photographs of the lean flame (ϕ = 0.8) without (left) and with (right) NRP discharges applied at 20 kHz. Figure extracted from [21]

  [START_REF] Blanchard | Experimental and numerical characterization of a lean turbulent premixed flame stabilized by nanosecond discharges[END_REF]

Fig. 3 :

 3 Fig. 3: Numerical solution of the 2D mean axial velocity on the centerline plane of the Mini-PAC configuration (case NP). The burner geometry and the positions of the electrodes are indicated. The dashed white line shows the position of radial profiles analyzed in section 4.

Fig. 4 :

 4 Fig. 4: Temporal evolution of the total thermal flame power for cases NP (black) and P (red).

Fig. 5 :

 5 Fig. 5: Temporal evolution of the flame surface for cases NP (black) and P (red). Dashed lines represent the flame surface resolved at the LES filter scale, while solid lines include also the subgrid flame wrinkling contribution.

Fig. 6 :

 6 Fig. 6: Instantaneous snapshots taken at the same time instant t= 20 ms of c * = 0.5 iso-surface colored by the heat release for cases NP (left) and P (right).

Fig. 7 :

 7 Fig. 7: Comparison of mean flame heat release rate (top) and mean CH 4 source term (bottom) radial profiles 4.5 mm above the bluff body for cases NP (dashed black line) and P (red line).

Fig. 8 :

 8 Fig. 8: Comparison of mean temperature at 4.5 mm above the bluff body for cases NP (dashed black line) and P (red line).

Fig. 9 :

 9 Fig. 9: Comparison of mean species mass fractions (H 2 O, CO 2 , O 2 , CO) at 4.5 mm above the bluff body for cases NP (dashed black line) and P (red line).

Fig. 10 :

 10 Fig. 10: Comparison of mean species mass fractions (OH, O, H 2 , H) at 4.5 mm above the bluff body for cases NP (dashed black line) and P (red line).

Fig. 11 :

 11 Fig. 11: chemical pathways in the discharge zone during (left) and after (right) a pulse.

Fig. 12 :

 12 Fig. 12: CH4 consumption pathways in the flame front without NRP discharges (dashed black line) and with NRP discharges (red line).

Acknowledgments

This work was supported by the ANR programs PASTEC (ANR-16-CE22-0005-01) and granted access to the HPC resources from IDRIS, TGCC and CINES under the allocation A0052B10253 made by GENCI (Grand Equipement National de Calcul Intensif) and HPC resources from the Mesocentre computing center of CentraleSupelec and Ecole Normale Superieure Paris-Saclay supported by CNRS and Région Ile-de-France ( http:// mesocentre.centralesupelec.fr/ ).We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France. Victorien Blanchard, Nicolas Minesi and Christophe Laux from EM2C lab are acknowledged for providing the experimental data.Vincent Moureau and Ghislain Lartigue from CORIA lab are acknowledged for providing the YALES2 code.