C Nicolò

; C Périer

; M Prague

; C Bellera

; G Macgrogan

; O Saut

; S Benzekry

S Benzekry

Machine learning and mechanistic modeling for prediction of metastatic relapse in early-stage breast cancer

Keywords: Breast cancer, mechanistic model, machine learning, metastatic relapse, survival analysis

Purpose: For patients with early-stage breast cancer, prediction of the risk of metastatic relapse is of crucial importance. Existing predictive models rely on agnostic survival analysis statistical tools (e.g. Cox regression). Here we define and evaluate the predictive ability of a mechanistic model for the time to metastatic relapse. Methods: The data consisted of 642 patients with 21 clinicopathological variables. A mechanistic model was developed on the basis of two intrinsic mechanisms of metastatic progression: growth (parameter 𝛼) and dissemination (parameter 𝜇). Population statistical distributions of the parameters were inferred using mixed-effects modeling. A random survival forest analysis was used to select a minimal set of 5 covariates with best predictive power. These were further considered to individually predict the model parameters, by using a backward selection approach. Predictive performances were compared to classical Cox regression and machine learning algorithms.

Results: The mechanistic model was able to accurately fit the data. Covariate analysis revealed statistically significant association of Ki67 expression with 𝛼 (p=0.001) and EGFR with 𝜇 (p=0.009). Achieving a c-index of 0.65 (0.60-0.71), the model had similar predictive performance as the random survival forest (c-index 0.66-0.69) and Cox regression (c-index 0.62 -0.67), as well as machine learning classification algorithms. Conclusion: By providing informative estimates of the invisible metastatic burden at the time of diagnosis and forward simulations of metastatic growth, the proposed model could be used as a personalized prediction tool of help for routine management of breast cancer patients.

Introduction

Breast cancer is the most frequent and second leading cause of cancer death in women 1 .

In the majority of cases, the disease is diagnosed at the early stage, when all detectable lesions, confined to the breast or nearby lymph nodes, can be surgically removed [START_REF] Noone | SEER cancer statistics review, 1975-2015[END_REF] .

However, approximately 20-30% of patients are reported to relapse with distant metastases after surgery [START_REF] Kohn | Invasion and Metastases[END_REF][START_REF] Pollard | Defining Metastatic Cell Latency[END_REF] , suggesting that clinically occult micro-metastases might already be present at the time of surgery. Accurate prediction of the risk of metastatic relapse is critical to personalize adjuvant treatment and avoid use of toxic and costly therapies when not needed.

In the era of artificial intelligence, prognostic models are playing an increasing role for such a task [START_REF] Shachar | Internet tools to enhance breast cancer care[END_REF] . Online tools, such as the Adjuvant! [START_REF] Ravdin | Computer Program to Assist in Making Decisions About Adjuvant Therapy for Women With Early Breast Cancer[END_REF][START_REF] Mook | Calibration and discriminatory accuracy of prognosis calculation for breast cancer with the online Adjuvant! program: a hospital-based retrospective cohort study[END_REF] and PREDICT models [START_REF] Wishart | PREDICT: a new UK prognostic model that predicts survival following surgery for invasive breast cancer[END_REF] , compute individualized survival probabilities based on multivariate statistical analysis and integration of clinical variables (age, tumor size, histological grade, hormone receptor status and nodal involvement) [START_REF] Shachar | Internet tools to enhance breast cancer care[END_REF] . These tools, however, are based on agnostic statistical models, such as Cox regression [START_REF] Wishart | PREDICT: a new UK prognostic model that predicts survival following surgery for invasive breast cancer[END_REF][START_REF] Wu | Personalized Prognostic Prediction Models for Breast Cancer Recurrence and Survival Incorporating Multidimensional Data[END_REF] . More recently, machine learning algorithms have started to be used [START_REF] Kourou | Machine learning applications in cancer prognosis and prediction[END_REF] .

Although traditionally designed for classification or regression tasks, adaptations to survival analysis include elastic net for Cox models [START_REF] Simon | Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent[END_REF] or the random survival forests algorithm [START_REF] Ishwaran | Random survival forests[END_REF] .

Deep learning has also recently been proposed for survival prediction from genomic data sets [START_REF] Yousefi | Predicting clinical outcomes from large scale cancer genomic profiles with deep survival models[END_REF] . Nevertheless, few studies have so far investigated machine learning for prediction of breast cancer survival or recurrence [START_REF] Yousefi | Predicting clinical outcomes from large scale cancer genomic profiles with deep survival models[END_REF][START_REF] Kim | Development of novel breast cancer recurrence prediction model using support vector machine[END_REF][START_REF] Delen | Predicting breast cancer survivability: a comparison of three data mining methods[END_REF] .

Mechanistic modeling approaches -where biological knowledge is used to build a simulation model -have been developed to describe metastatic dynamics [START_REF] Hanin | Uncovering the natural history of cancer from post-mortem cross-sectional diameters of hepatic metastases[END_REF][START_REF] Haeno | Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Optimum Treatment Strategies[END_REF][START_REF] Koscielny | A simulation model of the natural history of human breast cancer[END_REF][START_REF] Newton | A Stochastic Markov Chain Model to Describe Lung Cancer Growth and Metastasis[END_REF][START_REF] Retsky | Computer simulation of a breast cancer metastasis model[END_REF] . However, none of these mechanistic models has yet been implemented as a personalized predictive tool of metastatic relapse [START_REF] Shachar | Internet tools to enhance breast cancer care[END_REF] . In previous work, we evaluated a mechanistic model of metastatic development [START_REF] Iwata | A dynamical model for the growth and size distribution of multiple metastatic tumors[END_REF] using experimental data from orthosurgical mouse models of breast cancer [START_REF] Benzekry | Modeling spontaneous metastasis following surgery: an in vivo-in silico approach[END_REF] . The model was able to describe longitudinal growth of the total metastatic burden. For human data, the model could fit size-dependent probability of 20-years metastatic relapse in a historical dataset of 2,648 breast cancer patients [START_REF] Benzekry | Modeling spontaneous metastasis following surgery: an in vivo-in silico approach[END_REF][START_REF] Koscielny | Breast cancer: relationship between the size of the primary tumour and the probability of metastatic dissemination[END_REF] .

In the current work, we build on our descriptive model to propose an actionable tool for individualized predictions of the time to metastatic relapse (TTR), as well as reconstruction of the past natural history and prediction of future evolution of the disease. To train and validate the model, we relied on a dataset containing TTR and 21 clinical/pathological characteristics for 642 early-stage breast cancer patients. We first show a random survival forest analysis [START_REF] Ishwaran | Random survival forests[END_REF] , which allowed us to select a restricted number of predictors of interest. We then present the main novelty of this work: the calibration -using mixed-effects learning [START_REF] Lavielle | Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools[END_REF] of the mechanistic model. We illustrate the possible value of the mechanistic approach by performing predictive simulations of the entire cancer history of real patients, calibrated from data available at diagnosis only. Finally, we compare our results with predictive performances of classification machine learning algorithms for 5-years metastatic relapse.

Methods

Description of the data

The consisted of data of 642 women diagnosed with primary operable invasive breast carcinoma treated at the Bordeaux Bergonié institute between 1989 and 1993. This dataset has been comprehensively analyzed using standard statistical tools (Cox regression) [START_REF] De Mascarel | Comprehensive prognostic analysis in breast cancer integrating clinical, tumoral, micro-environmental and immunohistochemical criteria[END_REF] .

Patients in this analysis did not receive any adjuvant hormone or chemotherapy. The time to metastatic relapse (TTR) was defined as the time from the date of diagnosis to the date of distant relapse. Patients with no metastasis were censored at the date of last news or death [START_REF] Gourgou-Bourgade | Guidelines for time-to-event end point definitions in breast cancer trials: results of the DATECAN initiative (Definition for the Assessment of Time-to-event Endpoints in CANcer trials)[END_REF] . Clinical/pathological variables available in the dataset included age at diagnosis, menopausal status, histological grade, pathological T (tumor size) and N (axillary lymph node status) stages, pathological tumor size, histological type and number of metastatic lymph nodes. In addition tissue microarray array analysis was performed from the tumor samples. Area percentage of staining was determined using immunohistochemistry for estrogen receptor (ER), progesterone receptor (PR), HER2, Ki67, CK56, EGFR, VIM, CD24, CD44, ALDH1, BCL2, E-Cadherin and Trio [START_REF] De Mascarel | Comprehensive prognostic analysis in breast cancer integrating clinical, tumoral, micro-environmental and immunohistochemical criteria[END_REF] . Tumors were classified as HER2 positive if the HER2 immunostaining showed 3+ intensity or from individual review for 2+ scores [START_REF] De Mascarel | Comprehensive prognostic analysis in breast cancer integrating clinical, tumoral, micro-environmental and immunohistochemical criteria[END_REF] .

Missing covariate values were imputed before model building using the missForest imputation algorithm [START_REF] Stekhoven | MissForest-non-parametric missing value imputation for mixedtype data[END_REF] . The percentage of missing data was less than 5% for all variables (Figure S1). Using 100 trees per forest (predefined setting value of missForest), continuous and categorical covariates were imputed with a 4.4% and 7.1% error, respectively. Institutional review board approval was obtained for this retrospective study in accordance with national laws.

Random survival forests analysis (RSF)

The RSF algorithm is an extension of Breiman's random forest for the analysis of rightcensored time-to-event data [START_REF] Ishwaran | Random Survival Forests for R[END_REF] . We utilized the RSF implementation of the randomForestSRC R package [START_REF] Ishwaran | randomForestSRC: Random Forests for Survival, Regression, and Classification[END_REF] . All RSF models were fitted using 1000 trees, with the logrank splitting rule [START_REF] Ishwaran | Random Survival Forests for R[END_REF] . The optimal values of the tuning parameters (number of variables to be sampled at each split and minimum number of data points in a terminal node) were selected to maximize the concordance index calculated on the out-of-bag data.

Impact of covariates on the TTR was assessed using the forest-averaged minimal depth [START_REF] Ishwaran | High-Dimensional Variable Selection for Survival Data[END_REF] , which quantifies the predictive value of a covariate in a tree by its distance from the root node to the first node where it is used to split (smaller minimal depth values correspond to more predictive covariates).

Mechanistic model of metastatic dissemination and growth

The individual primary tumor (PT) kinetics in individual 𝑖 were described by the Gompertz model:

𝑉 % & (𝑡) = 𝑒 , -
. -0123 45 -6 [START_REF] Mook | Calibration and discriminatory accuracy of prognosis calculation for breast cancer with the online Adjuvant! program: a hospital-based retrospective cohort study[END_REF] ,

where 𝑉 % & (𝑡) is the number of cells within the PT at time 𝑡 and 𝛼 & and 𝑏 & are the Gompertzian growth parameters. Under this model, the tumor will approach a theoretical upper limit 𝐾 = 𝑒 ; -5 -. In order to avoid overfitting and improve identifiability, 𝐾 was fixed to 10 1> cells, leaving

𝛼 & as the only free parameter driving growth. This value was used following statistical estimations and biological considerations in breast cancer [START_REF] Koscielny | A simulation model of the natural history of human breast cancer[END_REF][START_REF] Spratt | Decelerating growth and human breast cancer[END_REF] . The PT size -reported as a diameter in the data -was converted into number of cells assuming spherical shape and the assumption 1mm @ = 10 A cells [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] . All model simulations were performed in number of cells.

Considering a dissemination rate from the PT given by 22

𝑑 & C𝑉 % & D = 𝜇 & 𝑉 % & ,
the total number of metastases at time 𝑡 is

𝑁 & (𝑡) = F 𝑑 & 0𝑉 % & (𝑠)7 𝑑𝑠 = F 𝜇 & 𝑉 % & (𝑠)𝑑𝑠. I J I J
The individual parameter 𝜇 & is the per day probability for a PT cell to disseminate and establish a distant metastatic colony. Each metastasis was assumed to start from the volume 𝑉 J of a single cell and to grow at the same rate than the PT:

𝑔 & (𝑣) = M𝛼 & -𝑏 & log M 𝑣 𝑉 J RR 𝑣.
The state of the metastatic process was described by a function 𝜌 & (𝑡, 𝑣) representing the distribution of metastatic tumors with size 𝑣 at time 𝑡. It is given by the solution of the following transport equation [START_REF] Iwata | A dynamical model for the growth and size distribution of multiple metastatic tumors[END_REF] :

𝜕 I 𝜌 & (𝑡, 𝑣) + 𝜕 V 0𝑔 & (𝑣)𝜌 & (𝑡, 𝑣)7 = 0, 𝑡 ∈ (0, +∞), 𝑣 ∈ (1, +∞),
endowed with the boundary and initial conditions:

𝑔 & (1)𝜌 & (𝑡, 𝑉 J ) = 𝑑 0𝑉 % (𝑡)7 , 𝑡 ∈ (0, +∞),

𝜌 & (0, 𝑣) = 0, 𝑣 ∈ (𝑉 J , +∞).

Mechanistic modeling of the time-to-relapse

To calibrate the metastatic model on TTR data, we defined the theoretical time to relapse as illustrated in Figure 1. More precisely, assuming a value 𝑉 V&Y as detection threshold, the time 𝜏 V&Y for a tumor to reach this size was given from the assumption of Gompertzian growth, i.e.

𝜏 V&Y = - 1 𝑏 & log [1 - 𝑏 & 𝛼 & log 𝑉 V&Y & \.
In an analogous way, the time from the first cancer cell to the detection of the primary tumor,

\.

A visibility threshold 𝑉 V&Y of 5 mm in diameter was assumed to be the detectability limit at imaging.

Since metastases of size larger than 𝑉 V&Y at time 𝑡 must have been emitted in the time interval C0, 𝑡 -𝜏 V&Y & D, the number of visible metastases at time 𝑡 can be obtained by

𝑁 V&Y & (𝑡) = 𝑁 & C𝑡 -𝜏 V&Y & D.
The theoretical TTR was then defined as

𝑇𝑇𝑅C𝑉 ]&^_ & ; 𝛼 & , 𝜇 & D = c inf g 𝑡 > 0 | 𝑁 V&Y & C𝑡 ]&^_ & + 𝑡D ≥ 1 k if 𝑁 & C𝑡 ]&^_ & D ≥ 1, +∞ otherwise.
That is, the time elapsed from diagnosis to the appearance of the first visible metastasis, if, according to the model, at least a metastasis was emitted before diagnosis; otherwise it was considered as infinite.

Calibration of the mechanistic model using mixed-effects learning

The A constant error model was assumed on the log-transformed data:

log 𝑇 𝑖 = log 0𝑇𝑇𝑅C𝑉 𝑑𝑖𝑎𝑔 𝑖 ; 𝜓 𝑖 D7 + 𝑒 𝑖 , (1)
where 𝑒 & is the residual error following a normal distribution with mean 0 and variance 𝜎 > .

The log-transformation was used to ensure positive values of the TTR variable. The individual parameters were assumed to be log-normally distributed and a linear covariate model was used:

log 𝛼 & = log 𝛼 %u% + 𝛽 , ⋅ 𝑥 , & + 𝜂 ,, & 𝜂 , & ∼ 𝒩(0, 𝜔 , > ) log 𝜇 & = log 𝜇 %u% + 𝛽 } ⋅ 𝑥 } & + 𝜂 } & , 𝜂 } & ∼ 𝒩C0, 𝜔 } > D (2)
where The contribution of individual 𝑖 to the likelihood is then given by [START_REF] Commenges | Dynamical Biostatistical Models[END_REF] :

ℒ & (𝜃) = F ℎC 𝑇 ‡ & | | 𝜓 & ; 𝜃 D • -𝑆C 𝑇 ‡ & | | 𝜓 & ; 𝜃 D𝑝C𝜓 & ; 𝜃D𝑑𝜓 & ,
where 𝜃 = C𝛼 %u% , 𝜇 %u% , 𝛽 ,, 𝛽 } , 𝜔 , , 𝜔 } , 𝜎D are the population parameters to be estimated and 𝑝C𝜓 & ; 𝜃D is the probability density function of the individual parameters, defined by (2).

The likelihood function (product of the individual likelihoods ℒ & (𝜃)) was maximized using the Stochastic Approximation of Expectation-Maximization (SAEM) algorithm [START_REF] Lavielle | Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools[END_REF] implemented in the R saemix package [START_REF] Comets | Parameter Estimation in Nonlinear Mixed Effect Models Using saemix , an R Implementation of the SAEM Algorithm[END_REF] . Initial values of the parameters were 𝛼 J = 0.002, 𝜇 J = 4 × 10 21> , 𝜔 , J = 15, 𝜔 } J = 15, 𝜎 J = 0.6. Standard errors of the estimated parameters were obtained using 100 bootstrap samples, and significance of covariates was assessed using the Wald test.

The mean of the conditional survival functions was estimated as

Classification machine learning algorithms for prediction of 5-year distant metastasis-free survival

Machine learning algorithms were trained using the scikit-learn python package [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] . These Each model was evaluated using 100 replicates of 10-fold cross-validations, with folds created using a stratified random sampling strategy (to preserve class balance). We applied a standard scaling for each run and each set by removing the training folds mean and scaling to unit variance, for the algorithms that require a homogeneous scale (SVM, kNN). The random forests algorithm was also seeded with a different random initialization value each time.

Evaluation of predictive performances

Results

Random survival forest multivariate analysis

For machine learning analysis of the TTR data, we used the RSF algorithm, which allows right-censored data [START_REF] Ishwaran | High-Dimensional Variable Selection for Survival Data[END_REF] . In addition to evaluating its predictive power, we used this algorithm to identify variables most predictive of TTR. Covariates were ordered on the basis of minimal depth (Figure 2A) and selected by running a nested analysis (Figure 2B). The crossvalidated c-index improved as the number of variables increased, reaching a maximum of 0.67 (95% CI, 0.66-0.69) with 5 variables. We selected these 5 top variables for the optimal model and future analysis in the mechanistic model below. Calibration plots for 2-, 5-and 10-year demonstrated a good agreement between the observed and model-predicted probabilities of metastatic relapse (Figure 2C), with moderate overprediction of relapse in the higher risk groups.

Partial dependence plots for the selected covariates indicated strong nonlinear relationships between covariates and relapse probabilities, with a non-monotonic behavior for age and tumor size (Figure S2). Confirming previous results [START_REF] Bellera | Variables with time-varying effects and the Cox model: Some statistical concepts illustrated with a prognostic factor study in breast cancer[END_REF] , these plots suggest nonvalidity of the proportional hazards assumption. Indeed, if such assumption would hold, the probability of no relapse at any time would be monotonous as a function of any covariate value.

Calibration and validation of the mechanistic model

To offer better insights on the mechanisms of relapse, we developed a mechanistic model of the TTR. We first evaluated the ability of this model to describe the TTR data without using covariates, except for pathological tumor size, which is a variable encoded in the structural model. We asked whether we could describe inter-individual variability of TTR by means of population statistical distributions of the parameters 𝜇 (dissemination) and 𝛼 (growth) of the model. Estimates of the population parameters were obtained using the SAEM algorithm and are reported in Table 1. Both fixed and random effects were identified with satisfactory precision (relative standard error < 37%). Figure 3A

Mechanistic covariate analysis and predictive power of the mathematical model

We next tested the covariates selected by the RSF analysis in the mechanistic model. We built the covariate model using a backward elimination procedure, starting with the full model with all the preselected covariates on both parameters 𝛼 and 𝜇. The final model included Ki67 and CD44 on 𝛼, and EGFR on the dissemination parameter 𝜇 (Table 1). The crossvalidated c-index for this model was 0.65 (95% CI, 0.60-0.71). Calibration plots for 2-, 5-and 10-year outcomes demonstrated good predictive accuracy of the model (Figure 3B).

Nevertheless, similarly as the RSF model, risk of relapse was overestimated in high risk groups. For the classification task of prediction of 5-and 10-year relapse, comparable results

were obtained between the mechanistic model, the RSF, classification machine learning algorithms and Cox regression (Table 2).

Predictive simulations of the mechanistic model

The previous results allow to calibrate the mechanistic model parameters from variables available at diagnosis, by using only the covariate part in equation ( 2) and neglecting the remaining unexplained variance. We used this to simulate the natural cancer history for a number of representative patients of our dataset. For each patient, the population-level parameters 𝛼 %u% , 𝜇 %u% , 𝛽 Ÿ&A ,, , 𝛽 OE¡¢¢,, and 𝛽 £¤¥¦,} were calibrated from an independent training set that did not contain this patient (coming from the cross-validation procedure).

Simulations were then performed using a discrete version of the metastatic model [START_REF] Benzekry | Modeling spontaneous metastasis following surgery: an in vivo-in silico approach[END_REF] S1). Distinct levels of Ki67 cause distinct growth kinetics. Moreover, unlike patient 224, the tumor of patient 358 expresses EGFR, which is associated with a higher metastatic potential. Thus although patient 358's tumor is much younger, the total number of (invisible) metastases at surgery is predicted to be larger (35 vs 21 metastases in patient 224). Model predictions are also informative in the case of individuals who were censored at the last follow-up. For instance, patient 70 (Figure S4) was censored at 17.7 years after diagnosis. This is consistent with our model, which predicts that this patient was disease-free after PT resection and would never have relapsed (TTR = +∞).

Comparison with machine learning classification algorithms and Cox regression

We tested the predictive power of machine learning classification algorithms. These cannot account for right-censored data. Thus, for this part we focused on prediction of 5-years metastatic relapse (yes or no). Best performances were achieved by the random forest and logistic regression models (Table 2 and Figure S5A-D). However, owing to the low event rate (9.25%), positive predictive value and and F1 scores were low (Table 2). To improve these metrics, a balanced, downsampled version of the dataset was constructed. This significantly improved PPV and F1, as well as model calibration (Figure S5E-F). We also compared our results to classical Cox regression survival analysis, which was found to exhibit similar predictive power (Figure S6, Tables 2 andS2).

Discussion

We propose a mechanistic model for prediction of metastatic relapse after surgical intervention in patients with early-stage breast cancer which, for the first time, is able to simulate the pre-and post-diagnosis history of the disease from data available at diagnosis only. Notably, with a c-index of 0.65 (95% CI, 0.60 -0.71), the mechanistic model achieved similar performances as the RSF algorithm (95% CI, 0.66-0.69) and Cox regression (95% CI, 0.62-0.67). These were also comparable with the actual standard prognosis model of breast-cancer-specific survival Adjuvant! (c-index 0.71) [START_REF] Mook | Calibration and discriminatory accuracy of prognosis calculation for breast cancer with the online Adjuvant! program: a hospital-based retrospective cohort study[END_REF] , used for risk classification in the MINDACT trial [START_REF] Cardoso | 70-Gene Signature as an Aid to Treatment Decisions in Early-Stage Breast Cancer[END_REF] . Others also reported a similar c-index of 0.67 for prediction of relapse in untreated patients, using Cox analysis [START_REF] Parker | Supervised risk predictor of breast cancer based on intrinsic subtypes[END_REF] . Advanced deep learning algorithms did not outperform this predictive power (c-index 0.68) for prediction of survival, even though integrating genomic data [START_REF] Yousefi | Predicting clinical outcomes from large scale cancer genomic profiles with deep survival models[END_REF] . A recent study that considered recurrence (either local, regional or distant) reported superior predictive power (AUC 0.81 for prediction of relapse at 5-years, versus 0.73 in our analysis), which might be explained by the much larger data set (15, 314 patients) and inclusion of epidemiological data not available in our analysis.

Grounded on the biology of the metastatic process, our model provides insights not achievable by statistical analysis alone. First, one of the most important and well known predictor of metastatic relapse -tumor size [START_REF] Koscielny | Breast cancer: relationship between the size of the primary tumour and the probability of metastatic dissemination[END_REF] -is directly incorporated into the model as an input parameter, used for calculation of the tumor age. Second, our model allows to test whether covariates are associated with growth and/or dissemination. We found that Ki67 was positively associated with the proliferation parameter 𝛼, which aligns with the definition of Ki67 as a proliferation marker [START_REF] Dowsett | Assessment of Ki67 in breast cancer: recommendations from the International Ki67 in Breast Cancer working group[END_REF] , as well as other results confirming its predictive power [START_REF] Yerushalmi | Ki67 in breast cancer: prognostic and predictive potential[END_REF] .

On the other hand, the basal marker EGFR was found to be a significantly positive covariate for the dissemination parameter 𝜇. EGFR has previously been demonstrated as being a prognostic factor of relapse [START_REF] Cheang | Basal-like breast cancer defined by five biomarkers has superior prognostic value than triple-negative phenotype[END_REF] , yet it is interesting to note that correlation appears here on 𝜇, which is consistent with the known fact that basal-like breast cancers are metastatically more aggressive [START_REF] Rakha | Basal-like breast cancer: a critical review[END_REF] . Interestingly, the median value of 𝜇 was consistent with the value estimated in a previous work using data of metastatic relapse probabilities from a cohort of breast cancer patients (𝜇 %u% = 2.26 × 10 21> cell -1 .day -1 here vs 𝜇 %u% = 7 × 10 21> cell -1 .day -1 in [START_REF] Benzekry | Modeling spontaneous metastasis following surgery: an in vivo-in silico approach[END_REF] ).

Our analysis also confirmed the prognostic value of age at diagnosis, with younger patients having a higher risk of relapse [START_REF] Lian | The Impact of Young Age for Prognosis by Subtype in Women with Early Breast Cancer[END_REF] . However, we found that after 60 years-old, the risk of relapse was increasing (Figure S2). This nonlinearity might explain why age did not appear as significant in either our mechanistic or Cox analysis. Although not significant at the 0.05 threshold, our results suggest prognosis value of CD44 (p = 0.083 for association with 𝛼).

CD44 is a cellular protein which is used as a marker of breast cancer stem cells [START_REF] Al-Hajj | Prospective identification of tumorigenic breast cancer cells[END_REF] and is associated with metastasis [START_REF] Li | Beyond tumorigenesis: cancer stem cells in metastasis[END_REF][START_REF] Mcfarlane | CD44 increases the efficiency of distant metastasis of breast cancer[END_REF] . Since CD44 is involved in cell adhesion and invasion [START_REF] Sheridan | CD44+/CD24-breast cancer cells exhibit enhanced invasive properties: an early step necessary for metastasis[END_REF] , it is surprising that it emerges as associated with 𝛼 and not 𝜇. This might indicate limitations of our model to detect complex biological processes since inference is only made indirectly and using several simplifying assumptions, such as a unimodal lognormal distribution of the parameters in the population. This might also be due to the fact that we considered here CD44 as a continuous variable, whereas a standard marker is the percentage of CD44 + /CD24 -cells [START_REF] Al-Hajj | Prospective identification of tumorigenic breast cancer cells[END_REF] . Similarly, only EGFR as a continuous variable was investigated and not the more classical EGFR + /CK5/6 + marker [START_REF] Cheang | Basal-like breast cancer defined by five biomarkers has superior prognostic value than triple-negative phenotype[END_REF] . Indeed, the aim of the current study was to establish the methodology of using our mechanistic model as a predictive tool, and we favored first keeping continuous variables. We plan to perform more detailed examination of such clinical variables in forthcoming work, as well as study the predictive power of the model in well-established subgroups such as node-negative patients or patients stratified according to the current molecular classification 50 .

A major advantage and clinical relevance of the mechanistic model over standard statistical or machine learning models, is that it can be used to perform patient-specific simulations allowing to assess the extent of invisible metastases at the time of diagnosis and to predict future growth of metastases. In turn, this might aid selecting patients that will most benefit from extended adjuvant therapy (or conversely, patients who would need only a limited number of cycles), by performing individualized simulations of the future course of the disease under competing therapeutic strategies. However, it will be first required to develop and validate models integrating the effect of systemic adjuvant therapies [START_REF] Benzekry | Modelling the impact of anticancer agents on metastatic spreading[END_REF] . We also believe that our methodology could be applied to other cancer types where similar concerns occur about the use of adjuvant therapy to avoid metastatic relapse (e.g. lung or kidney cancer).

In addition, the novel approach we propose to mechanistically model time-to-event data could be used to extract biologically relevant information from such data, which -although ubiquitous in clinical oncology -are almost exclusively analyzed using agnostic statistical tools.

Our model represents a first attempt of a mechanistic, individual-level, predictive model of metastatic relapse and might be improved in a number of ways. For instance, unexplained variability remained important despite the inclusion of covariates, suggesting that biomarkers other than those tested might improve model predictions. In this regard, genetic expression signatures have been shown to have higher predictive power compared to standard histological and clinical variables alone [START_REF] Van De Vijver | A gene-expression signature as a predictor of in breast cancer[END_REF] . Our mechanistic model could also be refined by higher order phenomena such as dormancy, which has been proposed to explain recurrence occurring after many years from surgery [START_REF] Karrison | Dormancy of Mammary Carcinoma After Mastectomy[END_REF][START_REF] Uhr | Controversies in clinical cancer dormancy[END_REF] . Finally, to be applied in clinical practice, the model should be further evaluated on external data sets. For comparison purpose between time-to-event and classification models, prediction metrics performed on the entire data set are reported. In bold is the best score achieved for a given metric.

Figure and table legends

Supplementary material

Minimal Depth

• • • • • • • • • • • • • • • • • • • • • 0

  mathematical model for individual TTR was thus defined as a function of 𝑉 ]&^_ & and a vector of structural individual parameters 𝜓 & = C𝛼 & , 𝜇 & D. Let 𝑇 & denote the TTR for patient 𝑖.

  𝑥 , & and 𝑥 } & are vectors of subject-specific covariates, which might be identical, or partially or completely different. The statistical model for the observations (1) implicitly defines for each individual 𝑖 the probability density function of 𝑇 & conditionally on the individual parameters 𝜓 & . The corresponding survival and hazard functions can be obtained by: 𝑆C 𝑡 | | 𝜓 & ; 𝑉 ]&^_ D = -𝑑 𝑑𝑡 log 𝑆C 𝑡 | | 𝜓 & D = 1 where Φ and 𝜑 are the cumulative distribution and probability density functions of the standard normal distribution, respectively. Let 𝐶 & denote the time to death or last follow-up for individual 𝑖. With right-censoring the observations are C𝑇 ‡ & , 𝛿 & D, where 𝑇 ‡ & = min(𝑇 & , 𝐶 & ) and 𝛿 & = 𝕝 gŠ -‹OE -k is the indicator variable.

(

  𝜓 1 , … , 𝜓 -) and (𝑉 ]&^_ 1 , … , 𝑉 ]&^_ -) drawn respectively from the estimated distribution of the individual parameters and the values of 𝑉 ]&^_ in the data.

1 >1(

 1 included logistic regression, support vector machines (SVM), random forests (RF), k-nearest neighbors (kNN) and gradient boosting. Values of the models hyperparameters were selected to optimize the area under the ROC curve, using a 3-fold cross-validated grid search. The kernel used for SVM was the radial basis function 𝐾(𝑥, 𝑥 š ) = 𝑒 2›‖•2•š‖ with 𝛾 = 21=number of features) and regularization parameter 𝐶 = 1. For the classification algorithms (and for them only), because these are not able to manage censored data, patients with data censored before 5 years were removed from the dataset, leaving a total of 594 patients. A balanced version of the dataset was constructed by keeping all 55 patients with 5-years metastatic relapse and by adding an equal number of patients randomly selected from those who experienced metastatic relapse later than 5 years from diagnosis.

  compares the model estimation of the TTR survival function to the empirical Kaplan-Meier estimate. Despite a slight overestimation of the metastatic risk for shorter times, the model was able to capture the shape of the Kaplan-Meier estimate. To further verify the agreement between model and data, we also compared model and Kaplan-Meier curves for different values of the tumor size at diagnosis (FigureS3). Although the model remained within the Kaplan-Meier confidence interval in all cases, we observed that model predictions tended to be less accurate as the tumor size increased, with overestimation of the relapse risk for large PT sizes.

For patients 224 and

  358 (Figure 4A-B), the model predicted a time from the first cell to detection (pre-surgical history) of 17.3 years and 3.61 years, respectively. In addition, the model allowed to predict the metastatic size distribution at diagnosis (Figure 4C-D). Patient 224 was predicted to have 21 invisible metastases at the time of diagnosis. The largest metastasis contained 2.32 • 10 4 cells and the smallest only 1 cell. The former was predicted to have initiated 12.9 years after the first cell of the PT, i.e. 4.4 years before surgery. The model predicts relapse 6.18 years after diagnosis while true relapse occurred at 4.88 years. At diagnosis, patient 358 was predicted to have a total of 35 invisible metastases. The largest metastasis contained 1.16 • 10 4 cells and the smallest 1.36 cells. According to this simulation, the first metastasis in this patient was emitted 2.6 years after the PT onset, i.e. 1.01 years before surgery. The model predicted relapse 1.97 years after diagnosis while true relapse occurred at 3.06 years. The differences in PT and metastatic dynamics for these two patients are due to the different values of the covariates (Table

Figure 1 :

 1 Figure 1: Scheme of the mechanistic model

Figure 2 :

 2 Figure 2: Survival random forest and covariate selection

Figure 3 :

 3 Figure 3: Fit of the mechanistic model

Figure 4 :

 4 Figure 4: Simulations of the cancer history for two patients

Figure S1 :Figure S2 :Figure S3 :

 S1S2S3 Figure S1: Percentage of missing values in each variable

Figure S6 :

 S6 Figure S6: Calibration plots for the Cox model with the eight covariates selected through

  Figure 1

Figure 3 +

 3 Figure 3

Figure S1 .

 S1 Figure S1. Percentage of missing values in each variable

  and are reported in Figures 4 and S4. Covariate values used for prediction and resulting inferred individual parameters 𝛼 § & and 𝜇̂& are reported in TableS1.

Table 1 :

 1 Parameter estimates. Mechanistic mixed-effects models for the time to metastatic relapse were fitted to the data using likelihood maximization. Parameters 𝛼 (growth) and 𝜇 (dissemination) were assumed to follow lognormal distributions with fixed effects (typical value) 𝛼 %u% and 𝜇 %u% , standard deviation of the random effects 𝜔 , and 𝜔 } and standard deviation of the error model 𝜎 (see methods). A first version was fitted without covariates (top) and a second with inclusion of dependence of the individual parameters on individual covariates (bottom). r.s.e. : relative standard error. p-value refers to a Wald test for statistical significance of the covariate on the parameter value.

Table 2 : Prediction metrics for classification of 5-and 10-years metastatic relapse.

 2

Table S1 :

 S1 Covariates and resulting inferred parameter values for individual predictions.

TTR = time-to-relapse

Table

S2

: Cox regression using the first five covariates selected by minimal depth with the random survival forest model. HR = hazard ratio. 95% CI = 95% confidence interval.

  + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

	Metastasis-free survival	+ ++ + + + + + + + + + + + + + + + + + ++ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.6 0.8 1.0
		0	5	10	15	20
			Time from diagnosis (years)

Table 1

 1 Parameter Estimate r.s.e. (%) p-value

		log α pop	-6.34	12.6
	Model without covariates	log µ pop σ ω α	-26.8 0.542 3.37	3.68 28.4 36.4
		ω µ	3.78	15.9
		log α pop	-9.01	10.8
		β Ki67,α	0.093	29.6	0.001
		β CD44,α	0.017	57.7	0.083
	Model with	log µ pop	-25.9	4.4
	covariates	β EGFR,µ	0.053	38.1	0.009
		σ	0.606	24
		ω α	2.75	22.1
		ω µ	3.03	20.5

Table 2

 2

		Algorithm	AUROC Accuracy Sensitivity Specificity PPV NPV	F1
		Mechanistic model	0.73	0.68	0.75	0.67	0.19 0.96 0.30
		Random survival forest	0.73	0.69	0.64	0.70	0.18	0.95	0.28
		Random forest	0.75	0.66	0.71	0.66	0.18	0.96	0.28
	5 years	Logistic regression k-nearest neighbor	0.75 0.62	0.83 0.91	0.42 0.02	0.87 1.00	0.24 0.41	0.94 0.91	0.31 0.05
		Gradient boosting	0.71	0.90	0.11	0.98	0.36	0.92	0.17
		Support vector machine	0.64	0.87	0.09	0.95	0.15	0.91	0.11
		Cox	0.71	0.72	0.66	0.73	0.20	0.95	0.31
		Mechanistic model	0.67	0.67	0.62	0.68	0.30 0.89 0.41
	10 years	Random survival forest	0.69	0.62	0.71	0.60	0.28	0.90 0.41
		Cox	0.65	0.65	0.61	0.65	0.28	0.88	0.39

&:𝑡 ]&^_ &

& D = 𝑃C 𝑇 & > 𝑡 | | 𝜓 & D = 1 -Φ ' log 𝑡 -log 0𝑇𝑇𝑅C𝑉 ]&^_ & ; 𝜓 & D7 𝜎 | | | | | 𝜓 & ƒ, ℎC 𝑡 | | 𝜓 & ; 𝑉 ]&^_ &

Acknowledgments

The authors thank Institut Bergonié for permission to use the data from the Breast Cancer Database for this research. They also thank the "Big bang" association for funding support, in memory of Marie-Christine Masini.

This study was achieved within the context of the Laboratory of Excellence TRAIL ANR-10-LABX-57.

All models were internally validated using 10-fold cross-validation. For classification tasks, the models were evaluated using receiver operating characteristic (ROC) curve analysis and standard performance metrics for binary classification algorithms [START_REF] Hossin | A review on evaluation metrics for data classification evaluations[END_REF] .

Discrimination of time-to-event models was quantified by Harrell's c-index [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF] . Calibration of each model was examined graphically, by comparing the model predicted probabilities against the observed event rates [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF] .

Time


Antoine Seguin

Yann Bertho

Philippe Gondret

Influence of confinement on granular penetration by impact

Keywords: numbers: 45.70.-n, 45.50.-j, 83.80.Fg

We study experimentally the influence of confinement on the penetration depth of impacting spheres into a granular medium contained in a finite cylindrical vessel. The presence of close lateral walls reduces the penetration depth, and the characteristic distance for these lateral wall effects is found to be of the order of one sphere diameter. The influence of the bottom wall is found to have a much shorter range.

I. INTRODUCTION

The motion of a solid sphere falling under its own weight into an unbounded viscous fluid has been known for more than one century to be characterized by the limiting Stokes velocity, and the influence of a confinement by close walls on the sphere velocity is also known for many years both theoretically [1,2,3] (see Ref. [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF] for a review) and experimentally [START_REF] Mackay | [END_REF]6]. This wall influence, due to the long range of hydrodynamic forces for low Reynolds number flows, results in a larger fluid friction force on the sphere. The velocity of a sedimenting sphere is thus known to be reduced near a wall either parallel or perpendicular to the sphere motion by a correction factor that is linear with d/l at low d/l (d/l < 1), where d is the sphere diameter and l is the distance from the wall. In all these cases, the well-known equations of classical fluid mechanics serve as a reliable support to understand and predict the corresponding motions with a fluid friction always proportional to velocity at small Reynolds numbers.

Considering now the motion of a falling sphere into granular matter, the existence of a nonzero solid friction force at vanishing velocity explains that the sphere driven by its own weight cannot reach a constant limit velocity but stops suddenly at a given depth. As the precise rheology of granular matter is far from being well understood and does not benefit from well-accepted equations despite recent important progress [7], the predictions for the stopping distance are far from easy. Many studies have been performed on the penetration depth δ of an impacting sphere into a granular medium. It emerges from those studies [8,9,10] that δ follows a power law of the form δ/d ∝ (ρ/ρ g ) β (H/d) α , where ρ g is the density of the grains, d and ρ are, respectively, the diameter and the density of the sphere, and H the total falling distance covered from release to rest. The first power exponent is β ≃ 0.5, whereas the second, α, varies between 0.3 and 0.5 depending mainly on the range of impact velocities and slightly on the packing fraction of the layer [8,9,10]. Indeed, as the packing varies from dense to loose random packing, the packing density does not vary so much, but the contact network and force network do vary significantly.

The knowledge and the expression of the forces exerted by the granular medium, and responsible for the sphere deceleration and stop, are still the subject of an intense debate [11,12], with different terms of frictional and collisional origins. Recently, the effect of a solid wall normal to a slow disk motion was shown to account for an exponential increase of the force on the disk in the close vicinity of the wall [13,14].

To our knowledge, no study deals in detail with the influence of close walls on the granular penetration, except the first results of Refs. [15,16]. In the present paper, after checking the usual "unbounded" case, we investigate experimentally the influence of the walls, either parallel or perpendicular to the sphere motion, for a sphere dropped onto a granular layer. In some situations, a coupling between the effects of parallel and normal walls may exist, but to understand their respective roles, each is studied separately in this paper.

II. EXPERIMENTAL SETUP

The experimental setup is sketched in Fig. 1. The penetration depth δ is investigated by dropping a spherical projectile of diameter d and density ρ onto a fine granu- lar material in a cylindrical container of diameter D and height b. The granular medium consists of glass beads (density ρ g ≃ 2.5 × 10 3 kg m -3 ) slightly polydisperse in size, with a diameter range of 300-400 µm. Before each drop, the granular medium is prepared by gently stirring the grains with a thin rod. The container is then overfilled and the surface leveled using a straight edge. We have checked that this preparation leads to reproducible results with only small variations. The grain size is much smaller than the falling sphere diameter d so that the granular medium can be considered as a continuum medium.

Different sphere materials and sizes have been used to bring out the influence of the sphere density ρ and of the sphere diameter d on the penetration depth δ. Steel spheres are initially maintained by a magnet at a distance h above the granular surface. Non-metallic spheres are held by creating locally a vacuum at the top of the sphere. Both apparatus allow one to drop the spheres without any initial velocity nor spinning motion. The sphere is released directly above the center of the container and falls along the container axis. The penetration depth δ is then measured with a thin probe by locating the top of the sphere with a precision better than 1 mm. The impact speed is tuned by varying the drop height h from 0 to 0.5 m; the corresponding velocity v at impact, given by v = √ 2gh where g is the gravitational acceleration, thus ranges from 0 to 3 m s -1 . In the following, the results will be analyzed using the total falling distance H = h + δ (Fig. 1).

III. THE UNBOUNDED CASE

We first present briefly our results obtained in the usual unconfined case: the container is large enough (D = 190 mm ≫ d) so that the sphere is not affected by the surrounding walls and also high enough (b = 300 mm ≫ δ) to avoid any bottom wall effect as we shall see in the next sections. The influence of the different parameters on the penetration depth is shown in Fig. 2 for five different sphere materials with a density ranging from 1.14 × 10 3 (polyamide) to 14.97 × 10 3 kg m -3 (tungsten carbide) and different sphere diameters d ranging from 5 to 40 mm. Note that each point corresponds to an average of about ten experiments, and error bars corresponding to the standard deviation are displayed only when larger than the symbol size. The penetration depth is larger for larger falling distance and larger projectile density. The data are well fitted by power laws of the form

δ d = A ρ ρ g β H d α , (1)
with A = 0.37±0.01, β = 0.61±0.02 and α = 0.40±0.04. The values of α and β are in close agreement with the previous studies already mentioned [9,10].

IV. INFLUENCE OF LATERAL CONFINEMENT

We are now interested in the influence of a lateral confinement on the penetration depth, using cylindrical containers of large height (b = 300 mm) but with smaller diameters D ranging from 128 down to 24 mm. The penetration depth δ is displayed in Fig. 3 as a function of the total falling distance H, for a steel projectile of diameter d = 19 mm falling in the different vessels. The two data sets corresponding to the two largest containers (D = 128 mm and D = 190 mm) coincide, suggesting that for large enough vessel diameters the influence of the surrounding walls becomes negligible and the medium can be considered as unbounded in the radial direction. For smaller diameters D, data still fall into straight lines in the log-log plot of Fig. 3 but with smaller slope values suggesting that Eq. ( 1) remains valid with smaller α values. The A and α values extracted from the fit of the data of Fig. 3 by Eq. ( 1) are reported in Figs. 4 and5.

The inset of Fig. 4 In this plot, data are well fitted by the exponential law

α = α ∞ 1 -exp - D -d λ α d , (2)
where α ∞ ≃ 0.4 corresponds to the "infinite" (unbounded) case and λ α ≃ 0.8 characterizes the range of wall effects. The characteristic distance λ α d of lateral wall effects is thus found here a little smaller than one sphere diameter, which is rather small. Hence, for D/d 5, the surrounding walls play a key part and prevent a deep penetration of the projectile. For D/d 5, the container has a vanishing influence and the sphere reaches a limiting penetration depth independent of D.

The prefactor A follows the same kind of decrease for decreasing container diameters D (see Fig. 5), but we believe that the relevant parameter for the lateral confinement is here D/d rather than (D/d -1), as A must vanish for vanishing D/d. The A variation is quite well fitted by the exponential law

A = A ∞ 1 -exp - D λ A d , (3)
where A ∞ ≃ 0.37 corresponds to the unbounded case and λ A ≃ 0.7 characterizes the range of wall effects. Note that the λ α and λ A values are about the same as these two parameters characterize the same wall effect. Besides, the power exponent β for the density ratio [Eq. ( 1)] does not depend on the lateral confinement but may be considered as constant, as data points for tungsten carbide (⋄) twice denser than steel fall on the same curve in both Figs. 4 and5.

The reduction of the penetration depth by a lateral confinement is clearly due to an enhanced blocking effect by the walls (less radial dilatancy). The precise mechanisms responsible for this reduction remain to be understood. For instance, the nonlinear "pressure" evolution inside the granular packing due to the well-known Janssen effect, and valid both in static and dynamic situations [17], fails in explaining an enhanced force on the sphere in the confined case.

V. INFLUENCE OF NORMAL CONFINEMENT

Let us now discuss the influence of a normal confinement by a close bottom wall. We investigate this effect by varying the thickness b of the granular layer contained in a vessel of large enough diameter D to avoid the lateral wall effects discussed before. As we are interested in the sphere penetration inside the granular material, we restrict however our study to thick enough layers with b > d. The normalized penetration depth The penetration depth remains unchanged by the presence of the bottom wall until the projectile approaches at a very short distance from the wall: δ follows the unbounded curve until the sphere impacts the bottom wall indicated by the dashed lines in Fig. 6. This very short range effect is consistent with recent measurements of the force on a flat plate approaching in a quasi-static way the solid bottom boundary of a granular sample, indicating that the penetration resistance increases exponentially near the boundary [13,14]. At a high enough falling height, the sphere impacts the wall (a characteristic shock sound can be heard) and bounces back leaving a crater below. For thin enough layers as already studied in Ref. [18], the opened crater does not become filled enough by in-wards avalanches so that no grains are at the center and the final rest position of the sphere corresponds to δ = b as the sphere comes to rest on the bottom "dry" wall [see (▽) in Fig. 6]. For thicker granular layers, inwards avalanches are large enough to bring back grains down to the crater center so that the sphere comes back at rest at δ < b [see ( ) in Fig. 6]. In that case, the larger the falling height, the longer the rebound time and thus the time for avalanche, so that the thicker the layer of back grains at the crater center: this explains the decreasing value of δ for an increasing falling height in this part of the curve. Data for b/d = 2.6 (•) would also display such a behavior for larger falling heights.

VI. CONCLUSION

We have measured the penetration depth δ of spheres of diameter d impacting a granular medium contained in a cylindrical vessel of diameter D. This penetration depth may be strongly affected by the presence of surrounding walls. The presence of a bottom wall perpendicular to the sphere motion has a very short range influence on δ as the sphere is only affected when the distance from the wall has reached a few grain diameters. On the other hand, the presence of lateral walls parallel to the sphere motion has a larger range influence with a typical characteristic distance λ ≃ d. For D/d 5, the walls have a strong influence and reduce significantly the penetration depth. For D/d 5, they have a vanishing influence. The two (bottom and lateral) wall effects have been studied independently but it could be a hard coupling between them for geometrically strongly hindered grain motions.

FIG. 1 :

 1 FIG. 1: (Color online) Sketch of the experimental setup and notations.

4 FIG. 2 :

 42 FIG. 2: Normalized penetration depth δ/d as a function of the normalized total falling distance H/d. Experimental data for steel spheres (ρ ≃ 7.8 × 10 3 kg m -3 ) of different diameters: (▽) d = 5 mm, (⋄) d = 10 mm, ( ) d = 19 mm, (•) d = 20 mm and (△) d = 40 mm. (-) Power law fit (δ/d) ∝ (H/d) α with α = 0.4. The shaded region corresponding to H < δ is by definition not allowed. Inset: (δ/d)(H/d) -0.4 as a function of the density ratio ρ/ρg for spheres of diameter d = 20 mm and of different materials: ( ) polyamide, (•) Teflon, ( ) glass, ( ) steel, and ( ) tungsten carbide. (--) Power law fit (δ/d)(H/d) -0.4 ∝ (ρ/ρg) β with β = 0.61.

  displays α as a function of D for the steel sphere of diameter d = 19 mm (•). A constant plateau value appears at large enough diameters D (D 100 mm) with the constant value α ≃ 0.4. At smaller D, α decreases significantly with D and seems to tend towards zero when D approaches the sphere diameter d = 19 mm, as no deep penetration (larger than d/2) would be possible for D < d. Data for different sphere diameters do not collapse in the α vs D plot, but they do collapse in the α vs (D/d -1) plot, as shown in Fig. 4.

FIG. 3 :FIG. 4 :

 34 FIG. 3: Normalized penetration depth δ/d for a steel sphere (d = 19 mm) as a function of the normalized total falling distance H/d for different cylindrical vessel diameters: ( ) D = 24 mm, (▽) D = 35 mm, ( ) D = 40 mm, ( ) D = 50 mm, ( ) D = 62 mm, (•) D = 80 mm, ( ) D = 128 mm, and (△) D = 190 mm.

FIG. 5 :

 5 FIG. 5: Prefactor A of Eq. (1) as a function of D/d with the same notations as in Fig. 4. (-) Best fit by Eq. (3) with A∞ ≃ 0.37 and λA ≃ 0.7. Inset: same data as a function of D.

FIG. 6 :

 6 FIG. 6: Normalized penetration depth δ/d of a steel sphere (d = 19 mm) as a function of the normalized total drop height H/d, for different heights b of the granular layer: (▽) b/d = 1.2, ( ) b/d = 1.9, (•) b/d = 2.6, ( ) b/d = 15.8. Each horizontal dashed line corresponds to each initial height b of the granular layer. (-) Power law fit (δ/d) ∝ (H/d) α with α = 0.4, for the unbounded case (b/d = 15.8).

ACKNOWLEDGMENTS

We wish to thank M. Salloum for preliminary experiments and J.-P. Hulin and M. Rabaud for a thoughtful reading of the manuscript.


Văn Chiến

Gérard C Bui

Gérard Henry

Edmond Duchamp

Ngoc Hoang

Quoc Minh

H Hoàn

Karol A Ngo

Penson

Henry Edmond

V C Bui

G H E Duchamp

V Hoang

Ngoc Minh

Q H Ngo

Hai Phong

Phan Dang Luu

Kien An

Viet Nam

K A Penson

A local Theory of Domains and its (Noncommutative) Symbolic Counterpart

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1 Introduction

This work is partly the continuation of [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF] where it has been established that the polylogarithms, indexed by the r-tuples (s 1 , . . . , s r ) ∈ C r , are well defined locally by

Li s 1 ,...,s r (z

) := ∑ n 1 >...>n r >0 z n 1 n s 1 1 . . . n s r r , for |z| < 1, (1)
could be extended, in case (s 1 , . . . , s r ) ∈ N r + , to some series, over the alphabet X = {x 0 , x 1 } generating the monoid X * with the neutral element 1 X * [START_REF] Berstel | Rational series and their languages[END_REF]. More precisely, (i) we start to consider [START_REF] Ngoc | Summations of Polylogarithms via Evaluation Transform[END_REF]

∀w = x s 1 -1 0 x 1 . . . x s r -1 0 x 1 ∈ X * x 1 , Li w = Li s 1 ,...,s r , (2)
(ii) then extend 2 Li • as the ⊔⊔ -morphism (C X , ⊔⊔, 1 X * ) -→ (C{Li} w∈X * , ×, 1) by adding Li x 0 (z) = log(z). This morphism is injective and satisfies [START_REF] Hoang Ngoc Minh | Aspects combinatoires des polylogarithmes et des sommes d'Euler-Zagier[END_REF] ∀S, T ∈ C X , Li S ⊔⊔ T = Li S Li T .

(iii) For the sake of symbolic calculations, it is important that, on the one hand, these series should belong to some "computable spaces" and, on the other hand, that the new domain (a) be closed by shuffle products and (b) that the Li • correspondence should preserve the shuffle identity [START_REF] Costermans | Some Results à l'Abel Obtained by Use of Techniques à la Hopf, Global Integrability of Field Theories and Applications[END_REF].

To this end a theory of global domains was presented in [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF]]. Here we focus on what happens in the neighbourhood of zero, therefore, the aim of this work is manyfold. Let us highligh the many facets of this matter.

(i) Propagate the extension to local Taylor expansions 3 as in [START_REF] Berstel | Rational series and their languages[END_REF] and the coefficients of their quotients by 1z, namely the harmonic sums, denoted H • and defined, for any w ∈ X * x 1 , as follows 4 [START_REF] Ngoc | Finite polyzêtas, Poly-Bernoulli numbers, identities of polyzêtas and noncommutative rational power series[END_REF] Li w (z)

1 -z = ∑ N≥0 H π X (w) (N)z N , (4)
by a suitable theory of local domains which assures to carry over the computation of these Taylor coefficients and preserves the stuffle indentity, again true for polynomials over the alphabet Y = {y n } n≥1 , i.e. ∀S,

T ∈ C Y , H S T = H S H T , (5)
2 This paper uses extensively shuffle and stuffle products (noted ⊔⊔ and respectively). For readers unfamiliar with these subjects their definitions are recalled at the end of this text, see paragraph 2.2. 3 Around zero. 4 Here, the conc-morphism π X : (C Y , conc, 1 Y * ) -→ (C X , conc, 1 X * ) is defined by π X (y n ) =

x n-1 0 x 1 and π Y its inverse on Im(π X ). See [6,5,13] for more details and a full definition of π Y .

meaning that H • : (C Y , , 1 Y * ) -→ (C{H w } w∈Y * , ×, 1), mapping any word w = y s 1 . . . y s r ∈ Y * to

H w = H s 1 ,...,s r = ∑ N≥n 1 >...>n r >0 1 n s 1 1 . . . n s r r , (6)
is a injective -morphism [START_REF] Ngoc | Finite polyzêtas, Poly-Bernoulli numbers, identities of polyzêtas and noncommutative rational power series[END_REF].

(ii) Extend these correspondences (i.e. Li • , H • ) to some series (over X and Y , respectively) in order to preserve the identity5 [START_REF] Ngoc | Finite polyzêtas, Poly-Bernoulli numbers, identities of polyzêtas and noncommutative rational power series[END_REF] Li π X (S) (z)

1 -z ⊙ Li π X (T ) (z) 1 -z = Li π X (S T ) (z) 1 -z . ( 7
)
true for polynomials S, T ∈ C Y .

(iii) Taking the definition of polyetas as in (1) at z = 1 or in ( 6) at +∞, one sees that, for any

s 1 > 1, Abel's theorem, one has ζ (s 1 , . . . , s r ) = lim z→1 Li s 1 ,...,s r (z) (8)
= lim N→+∞ H s 1 ,...,s r (N) (9) = ∑ n 1 >...>n r >0 1 n s 1 1 . . . n s r r . (10)
However, this theorem does not hold in the divergent cases. and we will recall some regularization process based on the computation of a -character with polynomial values and specialize it to obtain a character [START_REF] Costermans | Some Results à l'Abel Obtained by Use of Techniques à la Hopf, Global Integrability of Field Theories and Applications[END_REF][START_REF] Costermans | Noncommutative algebra, multiple harmonic sums and applications in discrete probability[END_REF] γ

• : (Q Y , , 1 Y * ) -→ (Z [γ], ×, 1), (11)
where Z := span Q {ζ (s 1 , . . . , s r )} r≥1,s 1 ≥2,s 2 ,...,s r ≥1 .

(iv) To this end, we use the explicit parametrization of the conc-characters obtained in [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF] and the fact that, under stuffle products, they form a group. We show the linear independence of the Kleene stars (z k y k ) * and show that γ • provides a group morphism between the group of conc-characters (endowed with ) and that of Taylor series g (with radius R = 1) such that g(0) = 1. This morphism maps each star y * k precisely to 1

Γ y k (1 + z) = exp -∑ n≥1 ζ (kn) (-z k ) n n , for k ≥ 2. ( 12
)
and y * 1 to the classical inverse Gamma:

Γ -1 y 1 (1 + z) = Γ -1 (1 + z). (13)
We will prove that all these "new" functions are entire and linearly independant.

To summarize, the present work concerns the whole project of extending H • over a stuffle subalgebra of rational power series on the alphabet Y , in particular the stars of letters and some explicit combinatorial consequences of this extension.

Domains and extensions

All starts with the (multiindexed) polylogarithm defined, for |z| < 1, by [START_REF] Berstel | Rational series and their languages[END_REF]. It is (multi-)indexed by a list (s 1 , . . . , s r ) ∈ N r ≥1 which can be reindexed by a word

x s 1 -1 0 x 1 . . . x s r -1 0 x 1 ∈ X * x 1 . From this, introducing two differential forms ω 0 (z) = z -1 dz and ω 1 (z) = (1 -z) -1 dz, (14)
we get an integral representation of the functions (1) as follows6 [START_REF] Ngoc | Summations of Polylogarithms via Evaluation Transform[END_REF] Li

w (z) =                  1 H (Ω) if w = 1 X * z 0 ω 1 (s) Li u (s) if w = x 1 u z 1 ω 0 (s) Li u (s) if w = x 0 u and |u| x 1 = 0, w ∈ x * 0 z 0 ω 0 (s) Li u (s) if w = x 0 u and |u| x 1 > 0, w / ∈ x * 0 , ( 15
)
where Ω is the simply connected domain C \ (] -∞, 0] ∪ [1, +∞[), over which we consider the algebra of analytic functions, H (Ω), with the neutral element 1 H (Ω) . This provides not only the analytic continuation of [START_REF] Berstel | Rational series and their languages[END_REF] to Ω but also extends the indexation to the whole alphabet X , allowing to study the complete generating series

L(z) = ∑ w∈X * Li w (z)w (16)
and show that it is the solution of the following first order noncommutative differential equation

     d(S) = (ω 0 (z)x 0 + ω 1 (z)x 1 )S, (NCDE) lim z∈Ω,z→0 S(z)e -x 0 log(z) = 1 H (Ω) X , asymptotic initinial condition, (17)
where, for any S ∈ H (Ω) X , for term by term derivation, one gets [START_REF] Drinfrl'd | -On quasitriangular quasi-Hopf algebra and a group closely connected with Gal Q/Q[END_REF]

d(S) = ∑ w∈X * d dz ( S | w )w. (18)
This differential system allows to show that L is a ⊔⊔ -character [START_REF] Hoang Ngoc Minh | Aspects combinatoires des polylogarithmes et des sommes d'Euler-Zagier[END_REF], i.e.

∀u, v ∈ X * , L | u ⊔⊔ v = L | u L | v and L | 1 X * = 1 H (Ω) . (19)
Note that, in what precedes, we used the pairing • | • between series and polynomials, classically defined by, for T ∈ C X and P ∈ C X7

T | P = ∑ w∈X * T | w P | w , (20)
where, when w is a word, S | w stands for the coefficient of w in S. With this at hand, we extend at once the indexation of Li from X * to C X by

Li P := ∑ w∈X * P | w Li w = ∑ n≥0 ∑ |w|=n P | w Li w . (21)
In [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF], it has been established that the polylogarithm, well defined locally by [START_REF] Berstel | Rational series and their languages[END_REF], could be extended to some series (with conditions) by the last part of formula (21) where the polynomial P is replaced by some series.

As was said previously, we focus here on what happens in the neighbourhood of zero. Therefore, the aim of this paragraph concerns the two first points of Section 1. which we summarize here (i) Propagate the extension to local Taylor expansions 8 of polylogarithms and the coefficients of their quotients by 1z, namely the harmonic sums, by a suitable theory of local domains.

(ii) Extend these correspondences (i.e. Li • , H • ) to some series in order to preserve the identity (7).

Polylogarithms: from global to local domains

The map Li • in general has been extended to a subdomain of C X , called Dom(Li • ) (see [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF]). It is the set of series

S = ∑ n≥0 S n ,
(z) = ∑ n≥1 z n n k = Li x k-1 0 x 1 (z) = L(z) | x k-1 0 x 1 but for t ≥ 0 (real), the series (tx 0 ) * x 1 belongs to Dom(Li • ) iff 0 ≤ t < 1. Dom(Li) C rat X A C X Above A = C X ⊔⊔ C rat X and C rat X
is the set of rational series [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF].

This definition has many merits 9 and can easily be adapted to arbitrary (open and connected) domains. But this definition, based on a global condition of a fixed domain Ω, does not provide a sufficiently clear interpretation of the stable symbolic computations around a point, in particular at z = 0. One needs to consider a sort of "symbolic local germ" worked out explicitely. Indeed, as the harmonic sums (or MZV) are the coefficients of the Taylor expansion at zero of the convergent polylogarithms divided by 1z, we only need to know locally these functions. In order to gain more indexing series and to describe the local situation at zero, we reshape and define a new domain of Li around zero to Dom loc (Li • ). The first step will be provided by the following theorem.

Theorem 2.2 Let S ∈ C X x 1 ⊕ C1 X * such that S = ∑ n≥0 [S] n where [S] n = ∑ w∈X * ,|w|=n S | w w, ([S] n are the homogeneous components of S), we suppose that 0 < R ≤ 1 and that ∑ n≥0 Li [S] n is unconditionally convergent (for the standard topology) within the open disk | z |< R. Remarking that 1 1 -z ∑ n≥0 Li [S] n (z) is unconditionally convergent in the same domain, we set 1 1 -z ∑ n≥0 Li [S] n (z) = ∑ N≥0 a N z N .
Then, for all N ≥ 0,

∑ n≥0 H π Y ([S] n ) (N) = a N .
Proof. Let us recall that, for any w ∈ X * , the function

(1 -z) -1 Li w (z) is analytic in the open disk |z| < R. Moreover, one has 1 1 -z Li w (z) = ∑ N≥0 H π Y (w) (N)z N . Since [S] n = ∑ w∈X * ,|w|=n S | w w and (1 -z) -1 ∑ n≥0 Li [S]
n absolutely converges (for the standard topology 10 ) within the open disk |z| < R, one obtains 1

1 -z ∑ n≥0 Li [S] n (z) = 1 1 -z ∑ n≥0 ∑ w∈X * ,|w|=n S | w w Li w (z) = ∑ n≥0 ∑ w∈X * ,|w|=n S | w w Li w (z) 1 -z = ∑ n≥0 ∑ w∈X * ,|w|=n S | w w ∑ N≥0 H π Y (w) (N)z N = ∑ N≥0 ∑ n≥0 ∑ w∈X * ,|w|=n S | w wH π Y (w) (N)z N = ∑ N≥0 H π Y ([S] n ) (N)z N .
This implies that, for any N ≥ 0,

a N = ∑ n≥0 H π Y ([S] n ) (N).

✷

We will need the following combinatorial Lemma 2.3 For a letter "a", one has

| (a + ) ⊔⊔ m | a n = m!S 2 (n, m) (23)
∑ n≥0 m!S 2 (n, m) x n n! = (e x -1) m . ( 24
)
Proof. (a + ) ⊔⊔ m is the specialization of

L m = a + 1 ⊔⊔ a + 2 ⊔⊔ . . . ⊔⊔ a + m
to a j → a (for all j = 1, 2 . . .m). The words of L m are in bijection with the surjections

[1 . . . n] → [1 . . .m], therefore the coefficient (a + ) ⊔⊔ m |a n is exactly the number of such surjections namely m!S 2 (n, m). A classical formula 11 says that ∑ n≥0 m!S 2 (n, m) x n n! = (e x -1) m . (25)

✷

To prepare the construction of the "symbolic local germ" around zero, let us set, in the same manner as in [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF],

Dom R (Li) := {S ∈ C X x 1 ⊕ C1 X * | ∑ n≥0 Li [S] n is unconditionally convergent in H (D <R )} (26)
and prove the following:

Proposition 2.4 With the notations as above, we have:

(i) The map ]0, 1] → C X given by R → Dom R (Li) is strictly decreasing (ii) Each Dom R (Li) is a shuffle subalgebra of C X .
Proof.

(i) It is straightforward that he map R -→ Dom R (Li) is decreasing. Set now, with

x + 1 = x 1 x * 1 = x * 1 -1, S(t) = ∑ m≥0 t m (x + 1 ) ⊔⊔ m
and let [S] n (t) be its homogeneous components, we have

∑ n≥0 Li [S] n (t) (z) = 1 -z 1 -(t + 1)z . For 0 < R 1 < R 2 ≤ 1 it is straightforward that Dom R 2 (Li) ⊂ Dom R 1 (Li).
Let us prove that the inclusion is strict. Take | z |< 1 and let us, be it finite or infinite, evaluate the sum

M(z) = ∑ n≥0 | Li [S] n (t) (z) |= ∑ n≥0 S(t) | x n 1 | Li x n 1 (z) |
then 11 See [START_REF] Stanley | Enumerative Combinatorics[END_REF], the twelvefold way, formula (1.94b)(pp. 74) for instance.

M(z) = ∑ n≥0 | S(t) | x n 1 | Li x n 1 (z) | = ∑ n≥0 ∑ m≥0 |t m (x + 1 ) ⊔⊔ m | x n 1 | Li x n 1 (z) | = ∑ m≥0 m!t m ∑ n≥0 S 2 (n, m) | Li x 1 (z) | n n! ≤ ∑ m≥0 m!t m ∑ n≥0 S 2 (n, m) Li n x 1 (| z |) n! , due to the fact that | Li x 1 (z) |≤ Li x 1 (| z |) (
Taylor series with positive coefficients). Finally, in view of equation ( 25), we get, on the one hand, for

| z |< (t + 1) -1 , M(z) ≤ ∑ m≥0 t m (e Li x 1 (|z|) -1) m = ∑ m≥0 t m ( | z | 1-| z | ) m = 1-| z | 1 -(t + 1) | z | .
This proves that, for all r ∈]0,

1 t + 1 [, ∑ n≥0 Li [S] n (t) (z) r < +∞.
On the other hand, if (t + 1) -1 ≤| z |< 1, one has M(|z|) = +∞, and the preceding calculation shows that, with t choosen such that

0 ≤ 1 R 2 -1 < t < 1 R 1 -1, we have S(t) ∈ Dom R 1 (Li) but S(t) / ∈ Dom R 2 (Li) whence, for 0 < R 1 < R 2 ≤ 1, Dom R 2 (Li) Dom R 1 (Li).
(ii) One has (proofs as in [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF]

) (a) 1 X * ∈ Dom R (Li) (because 1 X * ∈ C X ) and Li 1 X * = 1 H (Ω) . (b) Taking S, T ∈ Dom R (Li) we have, by absolute convergence, S ⊔⊔ T ∈ Dom R (Li). It is easily seen that S ⊔⊔ T ∈ C X x 1 ⊕ C1 X * and, moreover, that Li S Li T = Li S ⊔⊔ T 12 .

✷

In Theorem 2.6 bellow, we study, for series taken in C X x 1 ⊕ C.1 X * , the correspondence Li • to some H (D <R ), first (point 1) establishes its surjectivity (in a certain sense) and then (points 2 and 3) examine the relation between summability of the functions and that of their Taylor coefficients. For that, let us begin with a very general lemma on sequences of Taylor series which adapts, for our needs, the notion of normal families [START_REF] Montel | Leçons sur les familles normales de fonctions analytiques et leurs applications[END_REF]. Lemma 2.5 Let τ = (a n,N ) n,N≥0 be a double sequence of complex numbers. Setting

I(τ) := {r ∈]0, +∞[| ∑ n,N≥0 |a n,N r N | < +∞}, one has (i) I(τ) is an interval of ]0, +∞[, it is not empty iff there exists z 0 ∈ C \ {0} such that ∑ n,N≥0 |a n,N z N 0 | < +∞ (27)
In this case, we set R(τ) := sup(I(τ)), one has (a) For all N, the series ∑ n≥0 a n,N converges absolutely (in C). Let us note a Nwith one subscript -its limit (b) For all n, the convergence radius of the Taylor series Then I(τ) = / 0 and R(τ) ≥ R.

T n (z) = ∑ N≥0 a n,N z N is at least R(τ) and ∑ n∈N T n is summable for the standard topology of H (D <R(τ) ) with sum T (z) = ∑ n,N≥0 a N z N .
Proof.

(i) The fact that which proves that R(T n ) ≥ r, hence the result 13 . We also have

I(τ) ⊂]0, +∞[ is straightforward from the definition. If it exists z 0 ∈ C such that ∑ n,N≥0 | a n,N z N 0 |< +∞ then, for all r ∈]0, |z 0 |[, we have ∑ n,N≥0 | a n,N r N |= ∑ n,N≥0 | a n,N z N 0 | r | z 0 | N ≤ ∑ n,N≥0 |a n,N z N 0 | < +∞
| ∑ N≥0 a N r N |≤ ∑ N≥0 r N | ∑ n≥0 a n,N |≤ ∑ n,N≥0 | a n,N r N |< +∞ and this proves that R(T ) ≥ r, hence R(T ) ≥ R(τ).
(ii) Let 0 < r < r 1 < R and consider the path γ(t) = r 1 e 2iπt , we have

| a n,N |=| 1 2iπ γ T n (z) z N+1 dz |≤ 2π 2π r 1 T n K r N+1 1 ≤ T n K r N 1 with K = γ([0, 2π]), hence ∑ n,N≥0 | a n,N r N |≤ ∑ n,N≥0 | T n | K ( r r 1 ) N ≤ r 1 r 1 -r ∑ n≥0 T n K < +∞. ✷ Theorem 2.6 (i) Let T (z) = ∑ N≥0
a N z N be a Taylor series i.e. such that

lim sup N→+∞ | a N | 1/n = B < +∞, then the series S = ∑ N≥0 a N (-(-x 1 ) + ) ⊔⊔ N (28)
is summable (see [START_REF] Berstel | Noncommutative Rational Series with Applications, Encyclopedia of Mathematics and its Applications series[END_REF]) in C X (with sum in C x 1 ), S ∈ Dom R (Li) with R = (B + 1) -1 and Li S = T .

(ii) Let S ∈ Dom R (Li) and S = ∑ ≥0

[S] n (homogeneous decomposition), we define N -→ H π Y (S) (N) by 14 Li S (z)

1 -z = ∑ N≥0 H π Y (S) (N)z N . (iii) Moreover, ∀r ∈]0, R[, ∑ n,N≥0 | H π Y ([S] n ) (N)r N |< +∞, (29)
and, for all N ∈ N, the series (of complex numbers),

∑ n≥0 H π Y ([S] n ) (N) converges
absolutely to H π Y (S) (N). 13 For a Taylor series T , we note R(T ) the radius of convergence of T . 14 This definition is compatible with the old one when S is a polynomial.

(iv)

Conversely, let Q ∈ C Y with Q = ∑ n≥0
Q n (decomposition by weights), we suppose that it exists r ∈]0, 1] such that

∑ n,N≥0 | H Q n (N)r N |< +∞, (30)
in particular, for all N ∈ N,

∑ n≥0 H Q n (N) = ℓ(N) ∈ C unconditionally. Under such circumstances, π X (Q) ∈ Dom r (Li) and, for all z ∈ C, | z |≤ r, Li S (z) 1 -z = ∑ N≥0 ℓ(N)z N , (31)
Proof.

(i) The fact that the series (28) is summable comes from the fact that ω(a N (-(-x 1 ) + ) ⊔⊔ N ) ≥ N (see [START_REF] Berstel | Noncommutative Rational Series with Applications, Encyclopedia of Mathematics and its Applications series[END_REF]). Now from the lemma, we get

(S) n = ∑ N≥0 (a N (-(-x 1 ) + ) ⊔⊔ N ) n = (-1) N+n a N N!S 2 (n, N)x n 1 .
Then, with r = sup z∈K |z| (we have indeed r = Id K ) and taking into account that Li x 1 K ≤ log(1/(1r)), we have

∑ n≥0 Li (S) n K ≤ ∑ n≥0 ∑ N≥0 | a N | N!S 2 (n, N) Li x n 1 K ≤ ∑ n≥0 ∑ N≥0 | a N | N!S 2 (n, N) Li x 1 n K n! ≤ ∑ N≥0 | a N | ∑ n≥0 N!S 2 (n, N) | Li x 1 | n K n! ≤ ∑ N≥0 | a N | (e log( 1 1-r ) -1) N = ∑ N≥0 | a N | r 1 -r N .
Now if we suppose that r ≤ (B + 1) -1 , we have r(1r) -1 ≤ 1/B and this shows that the last sum is finite.

(ii) This point and next point are consequences of Lemma 2.5. Now, considering the homogeneous decomposition

S = ∑ n≥0 [S] n ∈ Dom R (Li).
we first establish inequation (29). Let 0 < r < r 1 < R and consider the path γ(t) = r 1 e 2iπt , we have

| H π Y ([S] n ) (N) |=| 1 2iπ γ Li [S] n (z) (1 -z)z N+1 dz |≤ 2π 2π Li [S] n K (1 -r 1 )r N+1 1 , K = γ([0, 1]
) being the circle of center 0 and radius r 1 . Taking into account that, for K ⊂ comp. D <R , we have a decomposition

∑ n∈N | Li [S] n | K = M < +∞, we get ∑ n,N≥0 | H π Y ([S] n ) (N)r N | = ∑ n,N≥0 | H π Y ([S] n ) (N)r N 1 | ( r r 1 ) N = ∑ N≥0 ( r r 1 ) N ∑ n≥0 | H π Y ([S] n ) (N)r N 1 | ≤ ∑ N≥0 ( r r 1 ) N M (1 -r 1 )r 1 ≤ M (1 -r 1 )(r 1 -r) < +∞.
The series

∑ n≥0 Li [S] n (z) converges to Li S (z) in H (D <R ) (D <R is the open disk defined by |z| < R).
For any N ≥ 0, by Cauchy's formula, one has,

H π Y (S) (N) = 1 2iπ γ Li S (z) (1 -z)z N+1 dz = 1 2iπ γ ∑ n≥0 Li [S] n (z) (1 -z)z N+1 dz = 1 2iπ ∑ n≥0 γ Li [S] n (z) (1 -z)z N+1 dz = ∑ n≥0 H π Y ([S] n ) (N)
the exchange of sum and integral being due to the compact convergence. The absolute convergence comes from the fact that the convergence of

∑ n≥ Li [S] n (z)
is unconditional [START_REF] Schaefer | Topological Vector Spaces[END_REF].

(iii) Fixing N ∈ N, from inequation (30), we get

∑ n≥0 | H Q n (N) |< +∞ which proves the absolute convergence. Remark now that (π X (Q)) n = π X (Q n ) and π Y (π X (Q n )) = Q n , one has, for all | z |≤ r | Li π X (Q n ) (z) |=| ∑ N∈N H Q n (N)z N |≤| ∑ N∈N H Q n (N)r N |, in other words Li π X (Q n ) D≤r ≤| ∑ N∈N H Q n (N)r N | and ∑ n∈N Li π X (Q n ) D≤r ≤| ∑ n,N∈N H Q n (N)r N |< +∞ which shows that π X (Q) ∈ Dom r (Li). The equation (31) is a consequence of point 2, taking S = π X (Q). ✷ Definition 2.7 We set Dom loc (Li) = 0<R≤1 Dom R (Li); Dom(H • ) = π Y (Dom loc (Li))
and, for S ∈ Dom loc (Li),

Li S (z) = ∑ n≥0 Li [S] n (z) and Li S (z) 1 -z = ∑ N≥0 H π Y (S) (N)z N .
Observe that, from this definition, theorem (2.8), will show that Dom(H

• ) is a stuffle subalgebra of C Y . (i) The series T = ∞ ∑ n=1 (-1) n-1 y n /n ∈ C Y is not in Dom(H • ) because, for all 0 < r < 1, one has ∑ n,N | T n (N)r N |≥ ∑ n≥0 1 1 -r = +∞ (32)
However one can get unconditional convergence using a sommation by pairs (odd + even).

(ii) For all s ∈]1, +∞[, the series

T (s) = ∞ ∑ n=1 (-1) n-1 y n n -s ∈ C Y is in Dom(H • ).
We can now state the Theorem 2.8 Let S, T ∈ Dom loc (Li), then

S ⊔⊔ T ∈ Dom loc (Li), π X (π Y (S) π Y (T )) ∈ Dom loc (Li)
and for all N ≥ 0,

Li S ⊔⊔ T = Li S Li T ; Li 1 X * = 1 H (Ω) , (33) H π Y (S) π Y (T ) (N) = H π Y (S) (N)H π Y (T ) (N). ( 34
) Li S (z) 1 -z ⊙ Li T (z) 1 -z = Li π X (π Y (S) π Y (T )) (z) 1 -z . ( 35
)
Proof. For equation (33), we get, from lemma 2.4 that Dom loc (Li) is the union of an increasing set of shuffle subalgebras of C X . It is therefore a shuffle subalgebra of the latter.

For equation (34), suppose S ∈ Dom R 1 0 (Li) (resp. T ∈ Dom R 2 0 (Li)). By [START_REF] Hadamard | Théorème sur les séries entières[END_REF] and theorem 2.6, one has

Li S (z) 1 -z ⊙ Li T (z) 1 -z ∈ Dom R 1 R 2 0 (Li),
where ⊙ stands for the Hadamard product [START_REF] Hadamard | Théorème sur les séries entières[END_REF]. Hence, for |z| < R 1 R 2 , one has

f (z) = Li S (z) 1 -z ⊙ Li T (z) 1 -z = ∑ N≥0 H π Y (S) (N)H π Y (T ) (N)z N (36)
and, due to theorem 2.6 point (iii), for all N, ∑ p≥0

H π Y (S p ) (N) = H π Y (S) (N) and ∑ q≥0 H π Y (T q ) (N) = H π Y (T ) (N) (absolute convergence) then,
H π Y (S) (N)H π Y (T ) (N) = ∑ p≥0 H π Y (S p ) (N) ∑ q≥0 H π Y (T q ) (N) = ∑ p,q≥0 H π Y (S p ) (N)H π Y (T q ) (N) = ∑ n≥0 ∑ p+q=n H π Y (S p ) π Y (T q ) (N) = ∑ n≥0 H (π Y (S) π Y (T )) n (N). (37)
Remains to prove that condition of Theorem 2.6, i.e. inequation (30) is fulfilled. To this end, we use the well-known fact that if ∑ m≥0 c m z m has radius of convergence R > 0, then

∑ m≥0 |c m | z m has the same radius of convergence (use 1/R = lim sup m≥1 | c m | -m ), then from the fact that S ∈ Dom R 1 0 (Li) (resp. T ∈ Dom R 2 0 (Li))
, we have (29) for each of them and, using the Hadamard product of these expressions, we get

∀r ∈]0, R 1 .R 2 [, ∑ p,q,N≥0 |H π Y (S p ) (N)H π Y (T q ) (N) r N | < +∞,
and this assures, for |z| < R 1 R 2 , the convergence of

f (z) = ∑ n,N≥0 H (π Y (S) π Y (T )) n (N)z N (38) applying theorem 2.6 point (iv) to Q = π Y (S) π Y (T ) (with any r < R 1 R 2 ), we get π X (Q) = π X (π Y (S) π Y (T )) ∈ Dom loc (Li) and f (z) = ∑ N≥0 ∑ n≥0 H (π Y (S) π Y (T )) n (N) z N = Li π X (π Y (S) π Y (T )) (z) 1 -z .
hence (34). ✷

Stuffle product and stuffle characters

For the some reader's convenience, we recall here the definitions of shuffle and stuffle products. As regards shuffle, the alphabet X is arbitrary and ⊔⊔ is defined by the following recursion (for a, b ∈ X and u, v ∈ X * )

u ⊔⊔ 1 X * = 1 X * ⊔⊔ u = u, ( 39
) au ⊔⊔ bv = a(u ⊔⊔ bv) + b(au ⊔⊔ v). ( 40
)
As regards stuffle, the alphabet is Y = Y N ≥1 = {y s } s∈N ≥1 and is defined by the following recursion

u 1 Y * = 1 Y * u = u, (41) y s u y t v = y s (u y t v) + y t (y s u v) + y s+t (u v). (42)
Be it for stuffle or shuffle, the noncommutative 15 polynomials equipped with this product form an associative commutative and unital algebra namely

(C X , ⊔⊔, 1 X * ) (resp. (C Y , , 1 Y * )).
Example 2.9 As examples of characters, we have already seen

• Li • from (Dom loc (Li • ), ⊔⊔ , 1 X * ) to H (Ω) • H • from (Dom(H • ), , 1 Y * ) to C N (arithmetic functions N -→ C)
In general, a character from a k-algebra16 (A , * 1 , 1 A ) with values in (B, * 2 , 1 B ) is none other than a morphism between the k-algebras A and a commutative algebra 17 B. The algebra (A , * 1 , 1 A ) does not have to be commutative for example characters of (C X , conc, 1 X * )i.e. conc-characters -where all proved to be of the form [START_REF] Duchamp | Towards a noncommutative Picard-Vessiot theory[END_REF]

∑ x∈X α x x * (43)
i.e. Kleene stars of the plane [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF]. They are closed under shuffle and stuffle and endowed with these laws, they form a group. Expressions like (43) (i.e. homogeneous series of degree 1) form a vector space noted C.Y . As a consequence, given P = ∑ i≥1 α i y i and Q = ∑ j≥1 β j y j , we know in advance that their stuffle is a conc-character i.e. of the form ( ∑ n≥1 c n y n ) * . Examining the effect of this stuffle on each letter (which suffices), we get the identity [START_REF] Duchamp | Towards a noncommutative Picard-Vessiot theory[END_REF]

∑ i≥1 α i y i * ∑ j≥1 β j y j * = ∑ i≥1 α i y i + ∑ j≥1 β j y j + ∑ i, j≥1 α i β j y i+ j * (44)
which suggests to take an auxiliary variable, say q, and code "the plane" C.Y , i.e. expressions like (43), like in Umbral calculus by

π Umbra Y : ∑ n≥1 α n q n -→ ∑ n≥1 α n y n (45)
which is linear and bijective 18 from

C + [[q]] to C.Y .
With this coding at hand and for S, T ∈ C + [[q]], identity (44) reads

(π Umbra Y (S)) * (π Umbra Y (T )) * = (π Umbra Y ((1 + S)(1 + T ) -1)) * (46)
This shows that if one sets, for z ∈ C and

T ∈ C + [[x]], G(z) = (π Umbra Y (e zT -1)) * (47)
we get a one-parameter stuffle group 19 , drawn on 1+C[z] + Y (a Magnus group), i.e. such that every coefficient is polynomial in z. Differentiating it we get

d dz (G(z)) = (π Umbra Y (T ))G(z) (48)
and (48) with the initial condition G(0) = 1 Y * integrates as

G(z) = exp (zπ Umbra Y (T )) (49)
where the exponential map for the stuffle product is defined, for any

P ∈ C Y such that P | 1 Y * = 0, is defined by exp (P) := 1 Y * + P 1! + P P 2! + . . . + P n n! + . . . . (50)
In particular, from (49), one gets, for k ≥ 1,

(zy k ) * = exp -∑ n≥1 y nk (-z) n n . (51)
This expression and that of

1 Γ(1 + z) = exp γz -∑ n≥2 ζ (n) (-z) n n ( 52
)
18 Its inverse will be naturally noted π Umbra q

. 19 i.e. G(z 1 + z 2 ) = G(z 1 ) G(z 2 ); G(0) = 1 Y * .
suggests to consider lacunary analogues of the inverse Gamma function together with a character which sends y 1 to γ and y n , n ≥ 2 to ζ (n). This -character is provided by asymptotic analysis of the Harmonic Sums. Indeed, one can show that, w ∈ Y * being given, the asymptotic expansion of N -→ H w (N), along the asymptotic scale (log(N) p N -q ) p,q∈N , at any rate 20 , can be written

∑ q≥0 Q w,q (log(N))N -q , ( 53
)
where Q w,q ∈ C[X ] (univariate polynomials) and, in particular,

Q w,0 ∈ Q[γ][X ] [3].
From this and the fact that H • is a -character, one gets that w -→ Q w,q (resp.

γ • : w -→ Q w,q (0)) is a -character with values in Q[γ][X ] (resp. Q[γ][X ]
). Now, a domain 21 Ω being given, it is easy to see that any -character χ (with general complex values and in particular γ • ) classically extends H (Ω) Y by

χ(P) = ∑ w∈Y * P | w χ | w (54)
as a -character from H (Ω) Y with values in H (Ω). Now, we can extend χ to some series, over Y . For that, let us set as above and as in [START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF], Definition 2.10 For any T ∈ H (Ω) Y , we note [T ] n the homogeneous compo-

nent 22 ∑ |w|=n T | w w of T Dom(χ, Ω) = {T ∈ H (Ω) Y |(χ(T n )) n∈N is summable in H (Ω)} (55)
The result, ∑ n≥0 χ(T n ) will be noted χ(T ).

This being defined, we have the following theorem ) 20 This means that the following expression is the limit of all partial asymptotic expansions. 21 Open, nonempty and connected subset of C. 22 The weight (w) of w ∈ Y * is just the sum of its indices 23 We will still note its extension to H (Ω) Y by χ. 24 In fact Dom(χ, Ω) is a subalgebra of (H (Ω) Y , ,

1 Y * )
Proof.

(i) By finitely supported sum.

(ii) S, T ∈ Dom(χ, Ω) then (χ(S p )) p≥0 , (χ(T q )) p≥0 are summable. But, as is graded for the weight, one has [S T ] n = ∑ p+q=n [S] p [T ] q . Take any K nonempty compact within Ω, then

∑ n≥0 χ([S T ] n ) K = ∑ n≥0 χ ∑ p+q=n [S] p [T ] q K = ∑ n≥0 ∑ p+q=n χ([S] p )χ([T ] q ) K ≤ ∑ n≥0 ∑ p+q=n χ([S] p ) K χ([T ] q ) K = ∑ p,q≥0 χ([S] p ) K χ([T ] q ) K < +∞.
The same computation without the seminorm proves (56).

(iii) If S ∈ Dom(χ, Ω) and we have to examine (and prove) the summability of the family (χ([exp (S)] n )) . Setting S = ∑ q≥0 [S] q , we have

[exp (S)] n = ∑ m≥0 ∑ q 1 +2q 2 +•••mq m =n [S] q 1 1 • • • [S] q m m q 1 !q 2 ! • • • q m ! . (57)
Hence, with all q i > 0,

∑ n≥0 χ([exp (S)] n ) K ≤ 1 + ∑ n>0 ∑ m>0 ∑ q 1 +2q 2 +•••mq m =n (58) χ([S] 1 q 1 K • • • χ([S] m ) q m K q 1 !q 2 ! • • • q m ! ≤ ∏ q≥1 e χ([S] q K = e ∑ q≥1 χ([S] q K = e M < +∞. (59)
because, as S ∈ Dom(χ, Ω), we have ∑ q≥1 χ([S] q ) K = M < +∞. ✷

A remarkable set of exponents

On the formal side, from (51), we have [START_REF] Duchamp | Towards a noncommutative Picard-Vessiot theory[END_REF] (z k y k ) * = exp

-∑ n≥1 y nk (-z) nk n , for z ∈ C, | z |< 1, (60)
as G r = G r otherwise it is an orbit as G r = ξ G r , where ξ is any solution of ξ r = -1 (this is equivalent to ξ ∈ G 2r and ξ / ∈ G r ). For r, q ≥ 1, we will need also a system X of representatives of G qr /G r , i.e. X ⊂ G qr such that

G qr = τ∈X τG r . (65)
It can also be assumed that 1 ∈ X as with X = {e 2ikπ/qr } 0≤k≤q-1 .

Proposition 2.13 ([5]) (i) For r ≥ 1, χ ∈ G r and z ∈ C, |z| < 1, the functions ℓ r and e ℓ r have the symmetry, ℓ r (z) = ℓ r (χz) and e ℓ r (z) = e ℓ r (χz) . In particular, for r even, as -1 ∈ G r , these functions are even.

(ii) For |z| < 1, we have (iii) For any odd r ≥ 2,

Γ -1 y r (1 + z) = e ℓ r (z) = Γ -1 (1 + z) ∏ χ∈G r {1}
e ℓ 1 (χz)

(iv) and, in general, for any odd or even r ≥ 2, ℓ r (z) = ∏ χ∈G r e ℓ 1 (χz) = ∏ n≥1

(1 + z r /n r ).

(v) For r ≥ 1, the function ℓ r is holomorphic on the open unit disc, D <1 , (vi) For r ≥ 1, the function e ℓ r (resp. e -ℓ r ) is entire (resp. meromorphic), and admits a countable set of isolated zeroes (resp. poles) on the complex plane which is expressed as χ∈G r χZ ≤-1 .

Proof. The results are known for r = 1 (i.e. for Γ -1 ). For r ≥ 2, we get (i) By (61), with χ ∈ G r , we get

ℓ r (χz) = -∑ n≥1 ζ (kr) (-χ r z r ) k k = -∑ k≥1 ζ (kr) (-z r ) k k = ℓ r (z),
thanks to the fact that, for any χ ∈ G r , one has χ r = 1. In particular, if r is even then ℓ r (z) = ℓ r (-z), i.e. ℓ r is even.

(ii) If r is odd, as G r = G r and, applying the symmetrization principle 26 , we get 26 Within the same disk of convergence as f , one has,

f (z) = ∑ n≥1 a n z n and ∑ χ∈G r f (χz) = r ∑ k≥1 a rk z rk . 22 -∑ χ∈G r ℓ 1 (χz) = -∑ χ∈G r ℓ 1 (χz) = r ∑ k≥1 ζ (kr) (-z) kr kr = ∑ k≥1 ζ (kr) (-z r ) k k .
The last term being due to the fact that, precisely, r is odd. If r is even, we have the orbit G r = ξ G r (still with ξ r = -1) and then, by the same principle,

-∑ χ∈G r ℓ 1 (χξ z) = r ∑ k≥1 ζ (kr) (-ξ z) kr kr = ∑ k≥1 ζ (kr) (-ξ z) r k k = ∑ k≥1 ζ (kr) -z r k k .
(iii) Straightforward.

(iv) Due to the fact that the external product is finite, we can distribute it on each factor and get

e ℓ r (z) = =1 ∏ χ∈G r e γ χz ∏ n≥1 χ∈G r 1 + χz n e -χz n = =1 ∏ n≥1 χ∈G r e -χz n ∏ n≥1 χ∈G r 1 + χz n .
Using the elementary symmetric functions of G r , we get the expected result.

(v) One has e ℓ 1 (z) = Γ -1 (1 + z) which proves the claim for r = 1. For r ≥ 2, note that 1 ≤ ζ (r) ≤ ζ (2) which implies that the radius of convergence of the exponent is 1 and means that ℓ r is holomorphic on the open unit disc. This proves the claim.

(vi) The function e ℓ r (z) = Γ -1 y r (1 + z) (resp. e -ℓ r (z) = Γ y r (1 + z)) is entire (resp. meromorphic) as finite product of entire (resp. meromorphic) functions, for r ≥ 1. The factorization in Proposition 2.13 yields the set of zeroes (resp. poles).

✷

As an example of projection, through γ• of an algebraic identity let us mention, for any z ∈ C, | z |< 1, one has, from (51), [START_REF] Duchamp | Towards a noncommutative Picard-Vessiot theory[END_REF] (z k y k ) * (-z k y k ) * = (-z 2k y 2k ) * , Γ y 2r (1 + z) = Γ y r (1 + ρz)Γ y r (1 + ρξ z),

where ρ is a 2r th -root of (-1) and ξ a primitive 2r th root of unity.

It is well known that the function e ℓ 1 (z) = Γ -1 (1 + z) is entire. In fact, all functions (64) are so (see Proposition 2.12). From this, ones get that (67) holds on the whole plane.

Example 2.14 [[12,11]] Let us give examples relating to polyzetas. For that, we use the following identities, for z ∈ C, | z |< 1, (see [START_REF] Duchamp | Towards a noncommutative Picard-Vessiot theory[END_REF]) (-zx 1 ) * ⊔⊔ (zx 1 ) * = 1 and (-zy 1 ) * ) * (zy 1 ) * = (-z 2 y 2 ) * , (-z 2 x 0 x 1 ) * ⊔⊔ (z 2 x 0 x 1 ) * = (-z 4 One can show (with a suitable extension of ζ , see [START_REF] Ngoc | Finite polyzêtas, Poly-Bernoulli numbers, identities of polyzêtas and noncommutative rational power series[END_REF][START_REF] Ngoc | Differential Galois groups and noncommutative generating series of polylogarithms[END_REF]) that 27 γ(-z 2 y 2 ) * = ζ ((-z 2 x 0 x 1 ) * ). Then, identifying the coeffients of z 2k , we get Again, with a suitable extension of 28 ζ (see [START_REF] Ngoc | Finite polyzêtas, Poly-Bernoulli numbers, identities of polyzêtas and noncommutative rational power series[END_REF][START_REF] Ngoc | Differential Galois groups and noncommutative generating series of polylogarithms[END_REF]) 27 Recall that, for any w ∈ Y * \ y 1 Y * , one has γ w = ζ (π X (w)) [START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF] and then it can be extended over series. 28 idem.

Conclusion

Noncommutative symbolic calculus allows to get identities easy to check and to implement. With some amount of complex and functional analysis, it is possible to bridge the gap between symbolic, functional and number theoretic worlds. This was the case already for polylogarithms and polyzetas. This is the project of this paper and will be pursued in forthcoming works.

(

  ii) Conversely, we suppose that it exists R > 0 such that (a) Each Taylor series T n (z) = ∑ N≥0 a n,N z N converges in H (D <R ). (b) The series ∑ n∈N T n converges unconditionnally in H (D <R ).

  in particular I(τ) = / 0 and it is an interval of ]0, +∞[ with lower bound zero. (a) Take r ∈ I(τ) (hence r = 0) and N ∈ N, then we get the expected result asr N ∑ n≥0 | a n,N |= ∑ n≥0 | a n,N r N |≤ ∑ n,N≥0| a n,N r N |< +∞. (b) Again, take any r ∈ I(τ) and n ∈ N, then ∑ N≥0 | a n,N r N |< +∞

Theorem 2 . 11

 211 Let χ : C Y -→ C be a -character 23 (i) H (Ω) Y ⊂ Dom(χ, Ω) (ii) If S, T ∈ Dom(χ, Ω) then S T ∈ Dom(χ, Ω) 24 and χ(S T ) = χ(S) χ(T )(56)(iii) If S ∈ Dom(χ, Ω), then exp (S) ∈ Dom(χ, Ω) and χ(exp (S)) = e χ(S

( 1 +

 1 ℓ r (z) = -∑ χ∈G r log(Γ(1 + χz)) and e ℓ r (z) = ∏ χ∈G r e γ χz ∏ n≥1 χz/n)e -χz/n .

  which, transformed by γ• and for | |z |< 1 and k ≥ 1, amounts to Euler's reflection formula (generalized to arbitrary k)[START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF]

x 2 0 x 2 1 )

 1 * and (-z 2 y 2 ) * (z 2 y 2 ) * = (-z 4 y 4 ) * .• From (67) and for r = 1, we haveγ(-z 2 y 2 ) * = γ(zy 1 ) * γ(-zy 1 ) * Γ -1 y 2 (1 + iz) = Γ -1 y 1 (1 + z)Γ -1 y 1 (1z) e -∑ n≥2 ζ (2n)z 2n /n = sin(zπ) zπ = ∑ k≥1 (ziπ) 2k (2k)! .

•

  Now, with r = 2, letting ρ 4 = -1, we have γ(-z 4 y 4) * = γ(z 2 y 2 ) * γ(-z 2 y 2 ) * Γ -1 y 4 (1 + z) = Γ -1 y 2 (1 + ρz)Γ -1 y 2 (1 + iρz),

γ(-z 4 y 4 )x 2 0 x 2 1 )

 41 * = ζ ((-z 4 y 4 ) * ), γ(-z 2 y 2 ) * = ζ ((-z 2 y 2 ) * ), γ(z 2 y 2 ) * = ζ ((z 2 y 2 ) * )then, using the poly-morphism ζ , we obtainζ ((-z 4 y 4 ) * ) = ζ ((-z 2 y 2 ) * )ζ ((z 2 y 2 ) * ) = ζ ((-z 2 x 0 x 1 ) * )ζ ((z 2 x 0 x 1 ) * )) = ζ ((-4z 4 * ).It follows then, by identification the coeffients of z 4k , that ζ (

  where S n := ∑

	Li k
		S | w	(22)
		|w|=n
	such that ∑ n≥0	Li S n is unconditionally convergent for the standard topology on H (Ω)
	[17].
	Example 2.1 [[10]] For example, the classical polylogarithms: dilogarithm Li 2 ,
	trilogarithm Li 3 , etc... are defined and obtained through this coding by

This research is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 101.04-2017.320

Here ⊙ stands for the Hadamard product[START_REF] Hadamard | Théorème sur les séries entières[END_REF].

Given a word w ∈ X * , we note |w| x 1 the number of occurrences of x 1 within w.

Here R is any commutative ring (like H (Ω), C, Z [γ], ...).

Around zero.

As the fact that, due to special properties of H (Ω) (it is a nuclear space[START_REF] Schaefer | Topological Vector Spaces[END_REF], see details in[START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF][START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF][START_REF] Ngoc | On the solutions of universal differential equation with three singularities[END_REF]), one can show that Dom(Li) is closed by shuffle products.

For this topology, unconditional and absolute convergence coincide[START_REF] Schaefer | Topological Vector Spaces[END_REF]

Proof by absolute convergence as in[START_REF] Duchamp | Kleene stars of the plane, polylogarithms and symmetries[END_REF].

For concatenation.

Here we will use k = Q or C.

In this context all algebras are associative and unital.

A family (g i ) i≥1 is said to be triangular if the valuation of g i , ϖ(g i ), equals i ≥ 1. It is easy to check that such a family is C-linearly free and that is also the case of families such that (g i -g(0)) i≥1 is triangular.

Acknowledgment

The authors would like to thank many colleagues for detailed reading of the manuscript and many valuable comments by which this paper could be improved.

and transform it through the -character γ• . First of all, we compute the radius of convergence of the image of the exponent (for coherence with the "bullet-notation", we will note γ y n the image of y n by the character γ • ) which gives [START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF] for

(61)

Then, from the fact that 1

are all in Dom( γ• , Ω) and therefore linearly independent. • and their transforms through γ• follow exponentiation (for z ∈ C, |z|< 1), i.e. [START_REF] Bui | Families of eulerian functions involved in regularization of divergent polyzetas[END_REF] γ(zy k

This leads us to set, for all k ≥ 1 and for

Proposition 2.12 ([5]) The families (ℓ r ) r≥1 and (e ℓ r ) r≥1 are C-linearly free and free from 1 H (Ω) .

Proof. Since (ℓ r ) r≥1 is triangular 25 then (ℓ r ) r≥1 is C-linearly free. So is (e ℓ re ℓ r (0) ) r≥1 , being triangular, we get that (e ℓ r ) r≥1 is C-linearly independent and free from 1 H (Ω) . ✷ Now, for any r ≥ 1, let G r (resp. G r ) denote the set (resp. group) of solutions, {ξ 0 , . . . , ξ r-1 }, of the equation z r = (-1) r-1 (resp. z r = 1). If r is odd, it is a group