Modibo Diabate
email: modibo.diabate@univ-grenoble-alpes.fr

Loren Coquille
email: loren.coquille@univ-grenoble-alpes.fr

Adeline Samson
email: adeline.leclercq-samson@univ-grenoble-alpes.fr

PARAMETER ESTIMATION AND TREATMENT OPTIMIZATION IN A STOCHASTIC MODEL FOR IMMUNOTHERAPY OF CANCER

Keywords: Immunotherapy, T cell exhaustion, Stochastic modeling, Mixed Effects Models, Treatment Optimization

Adoptive Cell Transfer therapy of cancer is currently in full development and mathematical modeling is playing a critical role in this area. We study a stochastic model developed by Baar et al. [4] for modeling immunotherapy against melanoma skin cancer. First, we estimate the parameters of the deterministic limit of the model based on biological data of tumor growth in mice. A Nonlinear Mixed Effects Model is estimated by the Stochastic Approximation Expectation Maximization algorithm. With the estimated parameters, we return to the stochastic model and calculate the probability of complete T cells exhaustion. We show that for some relevant parameter values, an early relapse is due to stochastic fluctuations (complete T cells exhaustion) with a non negligible probability. Then, focusing on the relapse related to the T cell exhaustion, we propose to optimize the treatment plan (treatment doses and restimulation times) by minimizing the T cell exhaustion probability in the parameter estimation ranges.

Introduction

Cancer is a group of more than 100 different diseases causing a large number of deaths a year worldwide [START_REF] Ferlay | Cancer incidence and mortality worldwide: sources, methods and major patterns in globocan 2012[END_REF]. Cells start growing uncontrollably due to genetic changes which impair their normal evolution. It can develop almost anywhere in the body. Cancer is a complex disease, difficult to study biologically (expensive and time consuming to experiment with animals and humans). In this context, mathematical modeling can be an excellent tool for stating or confirming biological assumptions with less expensive experiments.

When it is diagnosed quickly, cancer can be treated by chemotherapy, surgery, radiotherapy or by immunotherapy [START_REF] Airley | Cancer chemotherapy: basic science to the clinic[END_REF][START_REF] Hawley | Principles of radiotherapy[END_REF][START_REF] Liebelt | Principles of immunotherapy[END_REF][START_REF] Sabel | Principles of surgical therapy in oncology[END_REF]. Immunotherapy is a recent treatment that activates the immune system to kill cancer cells. In this paper we are interested by the Adoptive Cell Transfer (ACT) therapy to treat melanoma in mice with cytotoxic T cells, as experimented in 2012 by Landsberg et al. [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. This immunotherapeutic approach involves the stimulation of T cells which recognize one specific type of melanoma tumor cells (differentiated melanoma cells) through the antigen gp100 on their surface. The stimulated T cells are then able to kill these differentiated melanoma cells.

The authors of [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF] showed that during the inflammation induced by the therapy, proinflammatory cytokines called T N F α (Tumor Necrosis Factor) which are released in the body enhance a cell-type switch: the markers on the differentiated cancer cells disappear. They become dedifferentiated, and cannot get killed by T cells anymore. The resulting tumor tissue consists of both differentiated and dedifferentiated melanoma cells. Note that the switch is reversible (i.e. the melanoma cells can recover their initial type) and it does not require cell division or mutation. Figure 1 illustrates the described interactions between cells during the treatment.

In immunotherapy as in other cancer treatments, relapse is one of the main problems. The authors of [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF][START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF] describe two kinds of relapse in the ACT therapy. First, T cells only recognize the differentiated cancer cells (through the antigen gp100) and not the dedifferentiated cancer cells which they are thus not capable to kill. The growth of dedifferentiated cells induces a relapse. Second, T cells can become exhausted and no longer able to kill differentiated cancer cells thus causing a relapse. This problem was addressed in [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF] by T cells restimulation, which only lead to a delay in the occurrence of the relapse.

In these recent years, a lot of potentials have been seen in immunotherapy treatments, including the possibility of higher effectiveness [START_REF] Sagiv-Barfi | Eradication of spontaneous malignancy by local immunotherapy[END_REF] with lower side effects. But these promising new treatments are still to be understood, and long-term effects have not yet been studied in clinical trials, which are long and costly. An alternative is to study treatment failure and relapse through the development of mathematical models. Numerous deterministic models, mainly based on partial differential equations or ordinary differential equations, were developed for cancer study [START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF][START_REF] Novozhilov | Mathematical modeling of tumor therapy with oncolytic viruses: regimes with complete tumor elimination within the framework of deterministic models[END_REF][START_REF] Sun | A deterministic model of growth factor-induced angiogenesis[END_REF][START_REF] Watanabe | A mathematical model of tumor growth and its response to single irradiation[END_REF]. Stochastic counterparts include branching processes, individual-based or diffusion models [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF][START_REF] Cattiaux | A stochastic model for cytotoxic t. lymphocyte interaction with tumor nodules[END_REF][START_REF] Durrett | Branching Process Models of Cancer[END_REF][START_REF] Ochab-Marcinek | Pattern formation in a stochastic model of cancer growth[END_REF] and often a combination of both [START_REF] Anderson | Continuous and discrete mathematical models of tumor-induced angiogenesis[END_REF][START_REF] Quaranta | Mathematical modeling of cancer: the future of prognosis and treatment[END_REF]. The parameter estimation of those models often involve Mixed Effects Models [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS (statistics and computing)[END_REF] which are well adapted to observed longitudinal data [START_REF] Brant | Screening for prostate cancer by using random-effects models[END_REF][START_REF] Desmée | Nonlinear mixed-effect models for prostate-specific antigen kinetics and link with survival in the context of metastatic prostate cancer: a comparison by simulation of two-stage and joint approaches[END_REF][START_REF] Honerlaw | Biopsychosocial predictors of pain among women recovering from surgery for endometrial cancer[END_REF].

The purpose of this paper is to provide a quantitative study of the relapse and failure of skin cancer immunotherapy in mice using a stochastic individual-based model introduced in Baar et al. [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] and real biological data provided in Landsberg et al. [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. We first perform a statistical estimation of the parameters of the model, and then study the treatment relapse and failure focusing on T cell exhaustion. Finally we propose an optimization of the treatment plan.

The stochastic individual-based model recently proposed by [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] aims to describe the ACT therapy presented in the article [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. Its main features are: modeling the proliferation and the death of cancer cells and active T cells with birth and death processes [START_REF] Kendall | On the generalized" birth-and-death" process[END_REF] (the division of an active T cell is modeled as a birth event, and its exhaustion as a death event); modeling the switch between the two types of melanoma cancer cells; taking into account the interactions between cancer cells and T cells in a predator-prey framework (see related models in [START_REF] Berryman | The orgins and evolution of predator-prey theory[END_REF][START_REF] Costa | Stochastic eco-evolutionary model of a prey-predator community[END_REF]). The large population limit of the stochastic model is a deterministic differential system.

The authors of [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] have provided a set of biologically relevant parameters for which the stochastic system exhibits exhaustion of the T cells with high probability. However, these parameters are calibrated numerically but not estimated from real data. Moreover, they do not evaluate the quality of their model on real data, which is a key step if one wants to predict the tumor dynamics in a future step.

Our first objective is to validate the mathematical model, estimate its parameters and provide intervals of values by extracting the information from the biological data given in [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. We show that, in the range of parameters estimated from real data, the T cell exhaustion probability can indeed be non negligible. As the relapse caused by the exhaustion of T cells can be delayed by their restimulation, our second objective is to optimize the ACT protocol: we compute the optimal treatment doses and restimulation times by minimizing the T cell exhaustion probability at different stages of the disease evolution.

The paper is organized as follows. In Section 2 we present experimental data which we use for the parameter estimation. In Section 3, the stochastic model and its deterministic limit are described. In Section 4, we first present the statistical tools involved in the parameter estimation of the deterministic system, the method to estimate the T cell exhaustion probability and finally a procedure to optimize the ACT treatment. Section 5 presents the results: the estimated parameter values, the estimated exhaustion probability and the optimized treatment. The paper concludes with a discussion in Section 6.

Experimental data on Adoptive Cell Transfer therapy

We first describe the experimental setup of [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF], from which we use the data. Initially, melanoma cells are injected in each of the 19 mice such that the typical initial tumor diameter is less than 1 milimeter ( 10 5 cells). Mice are then split into 3 groups: untreated mice playing the role of a control group (denoted CTRL) and treated mice composed of two subgroups: mice treated once with ACT therapy (denoted ACT) and mice treated twice with ACT therapy (denoted ACT+Re). After a time t a = 70 days, mice of ACT and ACT+Re groups receive an intravenous delivery of 2 × 10 6 active T cells which are able to kill differentiated melanoma cells. We suppose that 2 × 10 3 stimulated T cells actually enter the tumor. The ACT+Re mice receive an additional injection of immunostimulatory nucleic acids in order to recover functions of exhausted T cells at time t Re = 160 days. In our modeling, both restimulated and injected T cells are considered as active T cells. Thus, the restimulation reactivates 2 × 10 3 T cells. Note that this value is not regarded as fixed; we study the impact of different stimulation and restimulation doses in Section 4.4.

During the experiment, the tumor development is measured by palpation when the tumor is small and digital photography when the tumor is larger. The tumor size is measured weekly using a vernier caliper and recorded as mean diameter. Longitudinal data for the 19 mice are composed of 5 mice of the CTRL group (with an average number of 10 observations per mouse); 7 mice of the ACT group (with an average number of 26 observations); and 7 mice of the ACT+Re group (with an average number of 33 observations).

The measure of tumor size, if it is smaller than 2 mm, is very inaccurate due to the difficulties of palpation. Thus, these values are left-censored. For tumor size between 2 and 3 mm, the inaccuracy is still important, about 1 mm. These values are intervalcensored in [2mm, 3mm]. Mice with tumor size exceeding 10 mm or showing signs of illness are killed. Thus measurements exceeding 10 mm are considered as right-censored. Tumor sizes between 3 and 10 mm are considered to be measured correctly with a measurement accuracy to be estimated.

Figure 2 represents the evolution of the tumor diameter along time for the 19 mice. The three censorships are illustrated by horizontal lines. Excluding the observations at t 0 = 0, mice in CTRL group are not concerned by left-censorship. Observe that the three groups of mice are distinguished by the speed of the tumor size evolution: tumor size of CTRL group (in blue) reaches 10 mm around the 100th day, tumor size of ACT group (in red) reaches 10 mm between 240 and 300 days and finally tumor size of ACT+Re group (in green) reaches 10 mm beyond the 300th day. Treatment and retreatment times are indicated by (black) vertical dotted lines. We distinguish two phases in the ACT therapy (for treated mice): the growth phase (before t a = 70) and the treatment phase (after t a = 70). Note that the mice in CTRL group only experience the tumor growth phase.

Modeling tumor growth under treatment

3.1.

A stochastic model with four cell types. The stochastic model proposed in [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] for the tumor growth under ACT therapy consists in a four-dimensional continuous time Markov process

Z(t) = (M (t), D(t), T (t), A(t)) ∈ N 4
where population sizes at time t are denoted by M (t) for differentiated melanoma cells, D(t) for dedifferentiated melanoma cells, T (t) for active T cells, and A(t) for cytokines T N F α . The carrying capacity K of the system is assumed to be equal to 10 5 cells, which corresponds to a tumor diameter of the order of the milimeter. This parameter scales the interaction rates between different cell populations, since the interaction is assumed to be mean-field (each cell interacts with all others, no space dependency is considered.) This value allows to study large impacts of stochastic events (such as complete exhaustion of T cells).

The main features of this stochastic individual-based process are the microscopic modeling of the proliferation of cancer cells and active T cells, the secretion of cytokines TNF-α, and the death, exhaustion, and clearance of, respectively, cancer cells, active T cells, and cytokines TNF-α. Those events are modeled with non-linear birth and death processes including a predator-prey interaction between cancer cells and T cells. The different events (division, death and switch events) occur at exponential random times whose rates constitute the parameters of the model. The complex biological process of T cell exhaustion is chosen to be modeled (and summarised) as the death of active T cells. More precisely, the time between the activation of a T cell and its exhaustion is assumed to be an exponentially distributed random variable.

Another mechanism which is not explicitly modeled concerns interferons IFN-γ which are a by-product of inflammation and strongly upregulate the expression of the markers on melanoma cells (see Supplementary Fig. 15 in [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]). Biologically, the upregulation of IFNγ-dependent markers and TNF-α-dependent dedifferentiation represent two functionally connected inflammation-induced adaptive mechanisms that together contribute to acquired resistance to ACT therapy [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. In the model of Baar et al. [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF], the role of interferons IFN-γ is qualitatively taken into account in the killing rate of differentiated melanoma cells by T cells. This significant simplification of real mechanisms underlying T cell exhaustion allows to work with a mathematically tractable model.

In the growth phase of the tumor (t < 70 days) only division, switch and death of melanoma cancer cells happen. At this stage, there are no active T cells nor cytokines T N F α : (T (t), A(t)) = (0, 0). The evolution of M (t) and D(t) is modeled by birth and death processes [START_REF] Kendall | On the generalized" birth-and-death" process[END_REF] including additional terms modeling the switches between the two cancer cell populations. Initial conditions are (M (0), D(0)) = (M 0 , 0).

At time t = 70 days (beginning of treatment), a dose of 2000 active T cells is injected in the system: T (70) = 2000. Then, all cell populations and cytokines (M (t), D(t), T (t), A(t)) evolve starting with initial conditions (M (70), D(70), 2000, 0). The system models the evolution of differentiated melanoma cells M (t) (prey) in the presence of active T cells T (t) (predators) as a predator-prey framework. An additional natural exhaustion of active T cells is modeled by a simple death process. A deterministic secretion of cytokines A(t) occurs at each T cell division event while the clearance of cytokines is modeled by a death process. An additional switch from differentiated to dedifferentiated melanoma cell due to the presence of cytokines is also modeled. The ACT group is modeled by the growth phase and this first treatment phase. At time t = 160 days, an additional dose of T-cell stimulant is injected for mice in ACT+Re group re-activating T cells. The quantity of reactivated T cells is unknown (contrary to the first injected dose) but it is biologically reasonable to assume that it activates the same quantity, i.e. 2000 T cells. Note that the impact of the two doses is studied in Section 4.4. Thus new initial conditions are (M (160), D(160), T (160) + 2000, A(160)). Figure 3 summarizes the dynamics of the process during the treatment phase (see Section 2.2 of [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] for more details). All events in the stochastic model (division, death and switch) occur after random exponential waiting times regulated by rate parameters. Parameters b M and b D (Figure 3) thus represent division rates of differentiated (respectively dedifferentiated) melanoma cells and correspond to the inverse of the average waiting time before observing the division of a differentiated (respectively a dedifferentiated) melanoma cell. Similarly, d M and d D correspond to natural death rates of differentiated (respectively dedifferentiated) melanoma cells, while s M D and s DM represent respectively the switching rates from differentiated melanoma cells to dedifferentiated one and the converse. 3 of the Appendix. 3.2. Deterministic limit of the stochastic model. The rescaled stochastic model converges in the limit of large population of cells (K → ∞) to the solution of the following deterministic differential system:

† A D M T d A d T s DM b D b M d D l prod A E b T M d M + d T M T s M D + s A A
       ṅM = (b M -d M )n M -s M D n M + s DM n D -s A n A n M -d M T n T n M ṅD = (b D -d D )n D + s M D n M -s DM n D + s A n A n M ṅT = -d T n T + b T n M n T ṅA = -d A n A + l prod A b T n M n T (3.1) with initial condition (n M 0 , n D 0 , n T 0 , n A 0 ). Quantity n X (t) ∈ R (with X = M, D, T or A)
represents the quantity of X at time t and ṅX represents the variation of this quantity along time. Note that we work in this paper on a simplified version of the initial model [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF]. Some competition terms were removed as they cannot be estimated from the data. Indeed, mice are killed before the tumor size (number of cancer cells) reaches the (non trivial) fixed point regulated by the competition terms. To avoid the problem of identifiability, we set

r M = b M -d M and r D = b D -d D since
only the difference between division and death rates of melanoma cells can be estimated from the data. The different phases of the therapy are similar to the ones of the stochastic model. During the growth phase (t < 70), only populations of cancer cells n M (t) and n D (t) evolve (n T (t), n A (t)) = (0, 0). All parameters associated with the treatment are then equal to zero. At time t = 70, the quantity of active T cells is set as n T (70) = D 70 . The four populations of cells (n M (t), n D (t), n T (t), n A (t)) evolve. Model (3.1) corresponds to the dynamics of the ACT group. For ACT+Re group, an additional quantity of active T cells is added in the system at time t Re = 160 (n T (160) = n T (t) + D Re ).

Methods

First, the statistical mixed model to estimate the parameters of the immunotherapy model from the data. A direct estimation of parameters is problematic: the likelihood function of the stochastic model is defined as a multiple integral over a high dimension event space. To the best of our knowledge, there exists no statistical method to estimate parameters of such models. As an approximation, we estimate parameters of the model through its deterministic limit using experimental data of tumor growth in mice provided by [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. We analyze these longitudinal data simultaneously with a Nonlinear Mixed Effects Model (NLMEM) [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS (statistics and computing)[END_REF] which takes into account the variability between mice, see Section 4.1. The parameters are estimated using a Stochastic Approximation Expectation Maximization (SAEM) algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF], see Section 4.2. Once the parameters have been estimated, T cell exhaustion is studied by returning to the stochastic model. Indeed, complete exhaustion is a purely stochastic phenomenon which is not modeled by the deterministic system. As this can be a rare event, we estimate its probability using a splitting Monte Carlo method [START_REF] Jacquemart-Tomi | A combined importance splitting and sampling algorithm for rare event estimation[END_REF][START_REF] Morio | An overview of importance splitting for rare event simulation[END_REF], see Section 4.3. Section 4.4 presents an optimization of the treatment in terms of treatment doses and restimulation times.

4.1. Mixed Effects Model for tumor growth under treatment. Parameters are estimated by analyzing the three groups simultaneously. We use a Nonlinear Mixed Effects regression Model (NLMEM) to take into account the inter-mice variability [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS (statistics and computing)[END_REF]. Parameter estimation with NLMEMs is difficult since the likelihood function does not have an explicit form. When the regression function of NLMEM is a stochastic model, the high dimension of the event space prevents to compute the likelihood. Thus, we choose to perform the parameter estimation based on the deterministic model (3.1).

Let us define y i = (y i1 , . . . , y in i ) where y ij is the (noisy) measurement of the tumor diameter for mouse i at time t ij , i = 1, . . . , N , j = 1, . . . , n i , and set y = (y 1 , . . . , y N ). The NLMEM is defined by:

y ij = f (ψ i , t ij ) + ij , (4.1)
ψ i = µ τ i exp{η i τ i },
where denotes component-wise multiplication and

µ = (r M , r D , s M D , s DM , n M 0 , b T , d T , d M T , s A , d A , l prod A ) is a vector of fixed effects; τ i = ( 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0 ) when i is a control individual; τ i = ( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) when i is a treated individual; , η i b T , η i d T , η i d M T , η i s A , η i d A , η i l A ) ∼ N (0, Ω)
η i is a vector of random effects independent of i ; ψ i represents the vector of individual regression parameters;

f (ψ i , t) = (n M (ψ i , t) + n D (ψ i , t)) 1 3
describes the tumor diameter (up to some constant), with n M and n D being solution of the nonlinear system (3.1); σ 2 is the residual variance; I n i the identity matrix of size n i ; Ω the variance matrix of the random effects quantifying variability between mice. We note θ = (µ, Ω, σ 2 ) the parameters to be estimated. Note that a biological constraint states that the natural switch rate from a dedifferentiated melanoma cell to a differentiated melanoma cell (s DM ) is larger than the rate in the other direction (s M D ) [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. This constraint is introduced in the estimation procedure. 4.2. Parameter estimation using SAEM algorithm. The SAEM-MCMC algorithm which combines the Stochastic Approximation Expectation Maximization algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] with a Markov chain Monte Carlo procedure adapted to censored data [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF][START_REF] Samson | Extension of the SAEM algorithm to leftcensored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF], implemented under the software MONOLIX [START_REF] Monolix | Monolix Users Guide Version 4.3.3, A software for the analysis of nonlinear mixed effects models[END_REF], is used to estimate the model parameters. We start with the model having random effects on all parameters and then use likelihood ratio tests (LRT) to select the significant random effects. Then, the effect of the categorical covariate group denoted by G is tested on the fixed effects µ through LRT. The ACT group is taken as reference group, G = {ACT*, ACT+Re, CTRL}. Standard errors of the estimated parameters of the selected models are estimated through the computation of the Fisher Information Matrix. 4.3. T cell exhaustion probability estimation using Monte-Carlo and Importance Splitting algorithms. T cell exhaustion is a cause of relapse [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF]. It is defined by the stochastic event {T (t) ≤ S, t ≤ t F } where S is the exhaustion threshold and t ≤ t F a condition for the exhaustion to occur in a finite time. The complete T cell exhaustion (S = 0) is a phenomenon which can only occur in the stochastic system: when the T cell population is low enough, the stochastic fluctuations can drive it to extinction, whereas the T cell population can never vanish in the deterministic limiting system. This follows from the analyticity of the solutions to (3.1), given that T (70) > 0.

We want to estimate the exhaustion probability p = P(T (t) ≤ S, t ≤ t F ). Let u = 1 T (t)≤S . Using the intuitive Monte Carlo method, the T cell exhaustion probability is estimated by p = N T k=1 u k /N T with N T the total number of simulations. However, depending on the values of the model parameters, the exhaustion threshold S may be difficult to reach leading to very small probabilities. Monte Carlo method thus requires a large N T when the probability p is very small. This method is therefore not suitable for estimating the probability of T cell exhaustion because the variance will diverge when the probability tends to zero for a reasonable value of N T [START_REF] Morio | An overview of importance splitting for rare event simulation[END_REF]. To reduce the variance for very small probabilities estimation, an alternative is the Importance Splitting (IS) algorithm designed for rare event probability estimation [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Jacquemart-Tomi | A combined importance splitting and sampling algorithm for rare event estimation[END_REF][START_REF] Morio | An overview of importance splitting for rare event simulation[END_REF]. Indeed, IS gradually calculates the probability of reaching the threshold S (the rare event) through the calculation of the probability of reaching intermediate thresholds easier to reach than S. Thus, the probability P(T (t) ≤ S, t ≤ t F ) is calculated according to the splitting principle by N T . To pass from iteration k to iteration k + 1, N T trajectories are sampled from the N S k trajectories having reached the threshold S k (by allowing replacement in the sampling) and run the N T new trajectories in order to reach the next threshold S k+1 . In our setting, T cells are completely exhausted at time t if the quantity T (t) has reached S 5 = 0. In Importance Splitting algorithms, it is possible to place the intermediate thresholds adaptively using a quantile method. We can also place these thresholds "manually" when we have enough knowledge about the dynamics of the stochastic system, which is the case for the dynamics of T cells. In this paper, we tested different numbers (3, 4 and 5) and spacing of intermediate thresholds until the numerical value of the exhaustion probability become stable from one series of 1000 simulations to another. We then fix these thresholds to S 1 = 7 × 10 -5 (or 0.35% of the initial quantity of active T cells), S 2 = 5 × 10 -5 (0.25% of T (70)), S 3 = 3 × 10 -5 (0.15% of T (70)), and S 4 = 2 × 10 -5 (0.1% of T (70)). We set N T to 1000.

P(T (t) ≤ S, t ≤ t F ) = m k=1 p k (4.2) with p 1 = P(T (t) ≤ S 1 , t ≤ t F ) and p k = P(T (t) ≤ S k | T (s) ≤ S k-1 , s ≤ t ≤ t F ) for k = 2, . . . ,
The computation of the T cell exhaustion probability requires N T = 1000 simulations of the stochastic model for each of the m = 5 thresholds. Furthermore, each trajectory can contain up to 100 millions of stochastic events. This stochastic approach is therefore very expensive in computing time. For example, the average time to calculate the exhaustion probability is approximately 100 hours (using 3 cores, with 7GB of memory per core).

However, this stochastic approach allows us to highlight a link between the probability of T cell exhaustion and the depth of the first minimum of the T cell deterministic trajectory. Indeed, since the stochastic process is locally close to the deterministic system plus a Gaussian noise [START_REF] Ethier | Markov processes: characterization and convergence[END_REF], the lower the deterministic minimum the larger the probability of T cell exhaustion. We exploit this relation to ease the study of the impact of remaining parameters which have a random effect on the T cell exhaustion phenomenon. 4.4. Criteria to optimize treatment doses and restimulation times. The kind of relapse we can act on without changing the medical setup is the relapse due to T cell exhaustion. The goal is thus to minimize the probability of T cell exhaustion, which amounts in the deterministic setting to maximize the value of the first minimum of the T cell population. Treatment parameters to be optimized are the treatment dose D 70 , the retreatment time t Re and the retreatment dose D Re . The cost function (to be maximized) is thus defined by

g = min{n T (t), t ≤ t F } (4.3)
with n T (t) the solution of deterministic system (3.1) and t F a fixed finite time. Note that g depends on all the model parameters through n T (t). Its unit is a number of T cells/K. We optimize the therapy by maximizing g over different subsets of treatment parameters. We refer to g as the minimum of the T cell trajectory, or simply, as the T cell minimum.

To take random effects into account in the optimization, we consider specific quantiles of their distribution. We focus on the T-cell related parameters d T (natural death rate of T cells), d M T (therapy induced death rate of the differentiated melanoma cells) and l prod A (which, multiplied by b T , leads to the therapy induced secretion rate of the cytokines). We then set the value of the natural clearance rate d A of cytokines to its population value since we observed numerically that the cost function g does not depend much on d A . We also set the value of the disease parameters r M and r D to their population mean assuming that these are disease-specific parameters that are less under control. As n M 0 corresponds to the initial size of the tumor, we consider different quantiles of its distribution. In summary, we set r M , r D and d A to their population values and consider quantiles of order 5%, 50%, 95% for d T , d M T , l prod A and n M 0 . These quantiles are denoted q α (φ i ), with α ∈ {5, 50, 95}% and the individual parameter φ i ∈ {d M T,i , d T,i , l prod A,i , n M 0 ,i }, and are provided in Table 1. Our choice of quantiles give a better idea on the dispersion of the estimated parameters in their estimation range compared, for example, to the quartiles. In total, we consider 3 4 = 81 sets of parameter values corresponding to three quantiles for each of the four parameters d M T , d T , l prod A , and n M 0 . We consider two optimization procedures: is sufficient to avoid T cell exhaustion before time t ≤ t F ). Indeed, from an experimental point of view, the initial dose d initial = 0.02 was chosen to have a good control on the tumor [START_REF]Private communication with Meri Rogava[END_REF]. This leads us to set D min , D max , D Re min and D Re max not too far from their initial values. Similar considerations lead us to set the retreatment time optimization interval to [t Re min , t Re max ] = [130,190]. Time t F is set to t F = 300 in the ACT group and to t F = 400 in the ACT+Re group (based on experimental data of Figure 2).

The optimal treatment parameter values (doses and restimulation time) are computed from the plots of the function g over intervals of discrete values of doses and restimulation time, for fixed sets of {d M T , d T , l prod A , n M 0 } taken in the estimation range.

Results

5.1.

No effect of the covariate group on the parameters. Likelihood ratio tests lead to the model given by {η

i r M = 0, η i r D = 0, η i l A = 0, η i d M T = 0, η i d T = 0, η i d A = 0, η i n M 0 = 0}.
LRT also showed that none of the parameters is function of the covariate group G, i.e., the model captures the differences between groups well enough without the help of an additional covariate. Thus, we take as final model the one without covariates. Table 1 contains the estimates and the associated standard errors of the final model, Figure 4 (and Figure 8 in Appendix) shows its individual fits. Figure 9 displays the observations versus predictions, Figure 10 the Visual Predictions Check (VPC), a prediction band obtained through simulations from the estimated model that is compared to real data [START_REF] Bergstrand | Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models[END_REF]. The VPC and the observations versus predictions plots are good, showing the good quality of the fit. The relative standard errors are small and in the range of what is expected for population parameter estimation [START_REF] Samson | The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model[END_REF], except for s M D . This might be an identifiability problem due to the indirect observation of the model, the four coordinates are not directly × 100 (in brackets). ω = 0 for parameters without random effect, and (-) when the standard error is not estimated. Quantile values for random effects considered in treatment optimization are given in the last two columns. Optimisation ranges of (re)stimulation doses and time are given at the bottom part of the table .   observed, only the third root of the sum of n M and n D . Eventually, confidence in the model predictions might be improved if model predictions could be related to experimentally measured parameters, like the expected ratios of T cells over tumor cells, or the ratios of antigen expressing over antigen loss tumor cells. However, we do not have access to such experimental data.

Note that the value of the death rate of T cells is very variable in different clinical cases (immunity against cancer, viruses, different markers involved, in-vitro or in-vivo clinical studies) because the processes underlying cell death are very different in each case. The values can range from a few days [START_REF] De Boer | Different Dynamics of CD4 + and CD8 + T Cell Responses During and After Acute Lymphocytic Choriomeningitis Virus Infection[END_REF] to a few weeks [START_REF] Pardoll | T cells take aim at cancer[END_REF] or months [START_REF] Westera | Closing the gap between T-cell life span estimates from stable isotope-labeling studies in mice and humans[END_REF] depending on the clinical situation. In our case, the half life of T cells is about 35 days, which is quite in accordance with [START_REF] Pardoll | T cells take aim at cancer[END_REF].

Moreover, the supplementary information of [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF] (Figure 9e) provides very partial information on the switch dynamics for in vitro experiments in the case of inhibited cell division and constant concentration of TNF-α. Those experimental data allowed Baar et al. to calibrate the switch parameters s M D , s DM , s A (see [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF], Supplementary Information, Section 2.3), and we checked a posteriori, using the same method as Baar et al. that our estimated parameter values are compatible with the data points. 5.2. Good fits for the three groups of mice. As shown in Figures 4 and8, data are well fitted, the growth and treatment phases are well enough captured despite the presence of censored data. The tumor growth phase is well fitted for CTRL mice as well as for treated mice (before the beginning of therapy). Treatment effect in the tumor size dynamics from the 70th day is also well highlighted in both ACT and ACT+Re groups. Furthermore, the restimulation effect is well marked in the tumor size dynamics for ACT+Re mice. 5.3. Biological relevance of the additional effect due to treatment s A . As cancer cells escape therapy, it is expected that the additional switching rate s A from differentiated to dedifferentiated melanoma cells is higher than the killing rate d M T of differentiated cells by T cells [START_REF]Private communication with Meri Rogava[END_REF]. Indeed, both mechanisms appear as quadratic terms involving the number of differentiated melanoma cells. Thus, if the cytokines and T cells populations sizes are similar, rates s A or d M T determine which mechanism prevails.

The estimated parameters confirm this relation: s A > d M T , as shown in Table 1. This is satisfactory concerning the biological modeling of the escape mechanism. Note that the calibration of model parameters in [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] led on the contrary to d M T much greater than s A . However, their results were calibrated manually without any comparison to real biological data. This reflects the relevance of our statistical approach. 5.4. Very small T cell exhaustion probability for population parameters. The T cell exhaustion probability estimated by the IS algorithm is very close to zero for individual and population parameters in both ACT and ACT+Re groups. Figure 5 shows the dynamics of cancer cells, T cells, and cytokines for population parameters in the ACT group. For these parameter values, T cells thus survive and continue to kill differentiated cells (whose number is controlled). However, the therapy enhances the switch towards dedifferentiated melanoma cells, which are not killed by T cells and thus grow exponentially. The tumor size regrows with a large ratio of dedifferentiated over differentiated melanoma cells. This is in accordance with Figure 1c of [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. 5.5. Evolution of T cell exhaustion probability with respect to therapy parameters. Some of the estimated parameters come with a non negligible random effect. We study the T cell exhaustion by naturally taking into account their distribution through the quantiles q α (φ i ).

Death rate of T cells d T is one of the most relevant therapy parameters. We thus estimate the T cell exhaustion probability as a function of d T in its estimation range (d T ∈ [0, q 99 (d T,i )]) when the other parameters are fixed to their estimated population values (still in ACT group). Figure 6 shows that the T cell exhaustion probability increases with d T , and even reaches non negligible values. This corroborates the conjecture in [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] asserting that there exist biological parameters for which the probability of relapse due to T cell exhaustion is non negligible. Note that when T cell exhaustion occurs, differentiated melanoma cells grow exponentially fast again, and become the most numerous in the final state of the tumor.

Figure 6 gives the probability of T cell exhaustion as a function of parameters with random effects (the others parameters being fixed to their population mean). We observe that the first minimum of the deterministic T cell trajectory becomes deeper and deeper as the parameters d T , d M T or l prod A increase. On the other hand, the minimum increases when disease parameters r M , r D or n M 0 increase. Furthermore, the minimum increases also with the cytokine clearance rate d A . Given these observations, we can expect that the joint variation of these parameters will lead us to interesting therapeutic scenarios. 5.6. Treatment optimization to minimize the probability of relapse due to T cell exhaustion. To reduce the number of parameter values of {d M T , d T , l prod A , n M 0 } to study, we focus on combinations of quantile values leading to an increased risk of T cell exhaustion. For that purpose, we define a threshold on the value of the function g (giving the value of the first minimum of the deterministic T cell trajectory) below which stochastic fluctuations have a good chance to lead to the T cell exhaustion. By setting this threshold to S 1 = 7 × 10 -5 (first threshold of IS algorithm), we identify 9 sets of parameter values of {d M T , d T , l prod A , n M 0 } for which the deterministic dynamics of the T cell population crosses the exhaustion threshold S 1 . Table 2 presents the optimization results for these 9 sets of parameter values presenting an increased risk of T cell exhaustion. Column ACT corresponds to optimization criterion (4.4) and column ACT+Re to criterion (4.5). Figures 7a and7b show optimal treatment parameters corresponding to two sets of parameter values for the optimization criterion (4.4) (Figure 7a presents a risk of T cell exhaustion, Figure 7b without an exhaustion risk). Figures 7c and7d show the optimization results for the same two sets of parameters but for criterion (4.5). Note that parameter values leading to an increased risk of T cell exhaustion mainly correspond to high values of parameters d T and l prod A as shown for their quantile of order 95% in Figure 7.

d M T , d T , l prod A , n M0 par. ACT g opt ACT g initial ACT ACT+Re g opt ACT+Re g initial ACT+Re q 5 , q 95 , q 95 , q 5 D opt
q 5 (d M T,i ), q 95 (d T,i ), q 50 (l prod A,i ), q 95 (n M 0,i ) Figure 7.
Global minimum of T cell population size according to treatment parameters using quantiles of order {5, 50, 95}%, for dMT , dT , l prod A , nM 0 . In ACT group (a,b): solid curve in cyan (q 5 (nM 0 ,i)), dashed curve in purple (q 50 (nM 0 ,i)), dotted curve in yellow (q 95 (nM 0 ,i)). Vertical lines: optimal doses D opt * 70 for the corresponding nM 0 . In ACT+Re group (c,d): curves for discrete retreatment times tRe = 130 (black solid), 142 (red dashed), 154 (green dotted), 166 (blue dotdashed), 178 (cyan longdashed), 190 (purple twodashed). Vertical lines: optimal D opt Re .

6. Discussion on optimized doses and restimulation times 6.1. Higher treatment doses and earlier restimulation times are not always better. Using quantiles of order 5%, 50%, 95% of d M T , d T , l prod A , and n M 0 , we observe different therapeutic scenarios. For most of parameter values, the optimal doses are the maximal ones: D opt 70 , D opt Re = D max (see Figure 7a). For other sets of quantile values (not presented in this paper), the optimal retreatment time is the minimal one t opt Re = t Re min . We can thus expect that higher treatment doses would further delay the T cell exhaustion or that an earlier restimulation would reduce the risk of an early exhaustion. However, the treatment parameters are optimized over optimization intervals which take into account biological considerations as explained in Section 4.4. Particularly, the Cytokine Release Syndrome (one of the most common, and potentially life threatening side effects of immunotherapy) is associated with increased doses of immunotherapy [START_REF] Maude | Managing cytokine release syndrome associated with novel t cell-engaging therapies[END_REF]. It may thus be experimentally non trivial to go beyond those intervals.

For some parameter values, the optimal restimulation dose is zero: D opt Re = 0. In these cases, the treatment dose D 70 is sufficient to avoid exhaustion. Finally, several parameter values (Figures 7b,7c,[START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF][START_REF] Desmée | Nonlinear mixed-effect models for prostate-specific antigen kinetics and link with survival in the context of metastatic prostate cancer: a comparison by simulation of two-stage and joint approaches[END_REF][START_REF] Durrett | Branching Process Models of Cancer[END_REF] led to optimal treatment doses and restimulation times which differ from the tested experimental ones. 6.2. Treatment optimization leads to values of g which are higher in the ACT+Re group than in the ACT group. Optimization benefit is evaluated by the probability of T cell exhaustion using the IS algorithm. As the function g is related to the exhaustion probability, we compare the values of g using respectively the initial treatment parameters (initial doses, and restimulation time provided in Section 4.4, and the optimized treatment parameters (Section 5.6). More precisely, let g initial G and g opt G (for G ∈ {ACT, ACT+Re}) represent the values of g G computed respectively using the initial treatment parameters and the optimized treatment parameters. Let Imp ACT+Re = (g opt ACT+Re -g initial ACT+Re )/g opt ACT+Re and Imp ACT = (g opt ACT -g initial ACT )/g opt ACT denote the relative shifts of g in the ACT+Re (resp. ACT) groups after treatment optimization, for the same values of the parameters. The two shifts are positive. Moreover, the difference Imp ACT+Re -Imp ACT is always positive with a mean value equal to (65.7 ± 4.8)% for the nine sets of parameters presenting an increased risk of T cell exhaustion (Table 2). Thus, the best improvements are noted in ACT+Re mice compared to ACT mice. This makes sense since ACT+Re mice have benefited from restimulation, which boosts the effects of optimization.

Note however that the treatment optimization only minimizes the risk of T cell exhaustion, the risk does not vanish completely as seen in Table 2: g opt is still smaller than S 1 . Considering wider optimization intervals could improve the results provided that the new limits remain biologically reasonable. 6.3. The larger n M 0 , the faster the tumor growth even for higher treatment doses. It has been biologically observed [START_REF]Private communication with Meri Rogava[END_REF] that when the tumor size exceeds a certain threshold at the beginning of the treatment (t = 70), the tumor control is less efficient even if more T cells are stimulated. = D max ). Large treatment doses D 70 do not slow down much the evolution of the tumor size for large initial value of differentiated melanoma cells. This is due to the fact that the initial ratio of tumour cells to T cells is higher when the tumours are larger : even the maximal dose does not allow to maintain this ratio constant.

Conclusion

In cancer treatment, resistance to therapy is a major problem. Understanding the mechanism of this resistance is then crucial. Parameter estimation for mathematical melanoma cancer immunotherapy model [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] using biological data from [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF] was essential in achieving this goal. Indeed, the estimated parameters lead to a statistically satisfactory model, which replaces the empirical calibration approach of Baar et al. [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF], and which constitutes a first necessary step for quantitative studies. We use these parameters to simulate realistic stochastic phenomena arising in the therapy and to highlight the existence of sets of biological parameters leading to each of the types of relapses identified by [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF]. Moreover, we are able to quantitatively estimate the probability of relapse due to T cell exhaustion, which allows to evaluate the quality of the treatment.

In addition, we confirme two important conjectures from the authors of [START_REF] Landsberg | Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation[END_REF]. First, the therapy escape mechanism suggests that the switch rate s A induced by cytokines should be higher than the rate d M T at which T cells kill differentiated melanoma cells. Second, biological observations [START_REF]Private communication with Meri Rogava[END_REF] show that when the tumor size exceeds a certain threshold at the beginning of the treatment, the tumor control is less efficient even if more T cells are stimulated. We indeed highlight the low effect of high treatment doses D 70 when the tumor size reaches a certain level at the beginning of the treatment.

With the estimated parameters, we optimize the treatment protocol using the deterministic system by plotting the value of the first local minimum of the T cell population as a function of treatment parameters. By linking the value of the minimum to the exhaustion probability of T cells, we provide optimal treatment doses and (re)stimulation time(s) which minimize the probability of relapse due to T cell exhaustion.

An alternative to the deterministic system would be a diffusion approximation of the stochastic model. A stochastic aspect would thus be included in the problem while keeping reasonable computational costs. An interesting statistical challenge, which is a work in progress, is to perform the parameter estimation in this stochastic diffusion approximation. The difficulties are mainly due to the predator-prey setting and the high dimensionality of the problem.

This paper focuses on the relapse mechanism of melanoma which is mediated by the T cell exhaustion. As explained in the introduction, we use a model which simplifies a lot the real mechanisms underlying the exhaustion process. Certain mechanisms must still be studied for a better understanding of the resistance to the therapy. This includes, for example, the role of IFN-γ, or the relapse due to dedifferentiated melanoma cells not killed by T cells. There are experimental prospects (requiring modifications of the treatment protocol) to handle the latter problem. One solution may be to stimulate another type of T cells which can kill the dedifferentiated melanoma (see [START_REF] Baar | A stochastic individual-based model for immunotherapy of cancer[END_REF] for a theoretical study of this treatment protocol). Another solution is to find a way to control the switching rate s A of differentiated melanoma to dedifferentiated melanoma cells, as recently considered by Glodde, Kraut et al. [START_REF] Glodde | Experimental and stochastic models of melanoma T-cell therapy define impact of subclone fitness on selection of antigen loss variants[END_REF]. Indeed, if s A is well controlled in the therapy (allowing the number of dedifferentiated cells to be driven to low enough values), then the problem of complete healing will come down to the control of T cell exhaustion. Thus, our work is a useful tool for new protocols of tumor treatment. T cell minimum according to treatment dose D70 in ACT group for dMT = q 5 (dMT,i). Top-figures: dT = q 5 (dT,i) ; Middle-figures: dT = q 50 (dT,i) ; Bottom-figures: dT = q 95 (dT,i). Left-figures: l prod A = q 5 (l prod A,i ) ; Middle-figures: l prod A = q 50 (l prod A,i ) ; Right-figures: l prod A = q 95 (l prod A,i ). In each sub-figure: nM 0 = q 5 (nM 0 ,i) (solid curve in cyan); nM 0 = q 50 (nM 0 ,i) (dashed curve in purple); nM 0 = q 95 (nM 0 ,i) (dotted curve in yellow). The vertical solid, dashed and dotted lines indicate the optimal dose D T cell minimum according to treatment dose D70 in ACT group for dMT = q 50 (dMT,i). Top-figures: dT = q 5 (dT,i) ; Middle-figures: dT = q 50 (dT,i) ; Bottom-figures: dT = q 95 (dT,i). Left-figures: l prod A = q 5 (l prod A,i ) ; Middle-figures: l prod A = q 50 (l prod A,i ) ; Right-figures: l prod A = q 95 (l prod A,i ).In each sub-figure: nM 0 = q 5 (nM 0 ,i) (solid curve in cyan); nM 0 = q 50 (nM 0 ,i) (dashed curve in purple); nM 0 = q 95 (nM 0 ,i) (dotted curve in yellow). The vertical solid, dashed and dotted lines indicate the optimal dose D . T cell minimum according to treatment dose D70 in ACT group for dMT = q 95 (dMT,i). Top-figures: dT = q 5 (dT,i) ; Middle-figures: dT = q 50 (dT,i) ; Bottom-figures: dT = q 95 (dT,i). Left-figures: l prod A = q 5 (l prod A,i ) ; Middle-figures: l prod A = q 50 (l prod A,i ) ; Right-figures: l prod A = q 95 (l prod A,i ).In each sub-figure: nM 0 = q 5 (nM 0 ,i) (solid curve in cyan); nM 0 = q 50 (nM 0 ,i) (dashed curve in purple); nM 0 = q 95 (nM 0 ,i) (dotted curve in yellow). The vertical solid, dashed and dotted lines indicate the optimal dose D T cell minimum according to treatment parameters in ACTRE group for dMT = q 5 (dMT,i); dT = q 95 (dT,i); l prod A = q 95 (l prod A,i ); nM 0 = q 5 (nM 0 ,i). From top to bottom and from left to right: D70 = 0.005, 0.012, 0.019, 0.026, 0.033, 0.040. In each subfigure: curves for discrete retreatment times tRe = 130 (black solid), 142 (red dashed), 154 (green dotted), 166 (blue dotdashed), 178 (cyan longdashed), 190 (purple twodashed). Vertical lines: optimal D opt Re . D T cell minimum according to treatment parameters in ACTRE group for dMT = q 5 (dMT,i); dT = q 95 (dT,i); l prod A = q 95 (l prod A,i ); nM 0 = q 95 (nM 0 ,i). From top to bottom and from left to right: D70 = 0.005, 0.012, 0.019, 0.026, 0.033, 0.040. In each subfigure: curves for discrete retreatment times tRe = 130 (black solid), 142 (red dashed), 154 (green dotted), 166 (blue dotdashed), 178 (cyan longdashed), 190 (purple twodashed). Vertical lines: optimal D opt Re . D Tumor size along time using quantiles q α (φi), α ∈ {5, 50, 95}%, φ ∈ {dT M , dT , l prod A , nM 0 } with D70 = Dmax, dT = q 95 (dT,i), l prod A = q 95 (l prod A,i ). Top: dMT = q 5 (dMT,i); Middle: dMT = q 50 (dMT,i); Bottom: dMT = q 95 (dMT,i). For each subfigure, nM 0 = q 5 (nM 0 ,i) in black solid line, nM 0 = q 50 (nM 0 ,i) in red dashed line, nM 0 = q 95 (nM 0 ,i) in green dotted line. Vertical orange line: treatment start.

Figure 1 .

 1 Figure 1. Cell interactions in ACT therapy: T cells kill differentiated melanoma cells and produce T cells and cytokines T N Fα (line 1); switch between differentiated and dedifferentiated melanoma cells is reversible (line 2); T cells cannot kill dedifferentiated melanoma cells (line 3).

Figure 2 .

 2 Figure 2. Tumor diameter (mm) along time (days) in the groups CTRL (dashed line in blue), ACT (solid line in red), ACT+Re (dotted line in green). Censorship indicated by horizontal lines. Vertical dotted lines: Treatment (day 70) and re-treatement (day 160) (for ACT+Re mice).

Figure 3 .

 3 Figure 3. Melanoma stochastic process: M : differentiated melanoma cells; D: dedifferentiated melanoma cells; T : T cells; A: cytokines T N Fα. †: death of cancer cells and clearance of cytokines; E: T cells exhaustion. Arrow with same initial and end points indicates a division in this population. Values above the arrows represent the rate parameters.

  Parameters b T and d T are respectively the division and exhaustion rates of active T cells. To a certain extent, b T can be seen as an effective parameter corresponding to b T = b 0 T -d 0 T with b 0 T and d 0 T respectively the therapy induced division and exhaustion rates of active T cells. Indeed, differentiated melanoma cells secrete a substance that exhausts active T cells with a rate d 0 T in addition to their natural exhaustion rate d T . The therapy induced death of differentiated melanoma cells is regulated by the rate parameter d M T while parameter d A corresponds to the clearance rate of cytokines T N F α . The product l prod A × b T represents the secretion rate of cytokines induced by the therapy, and s A is associated to the additional switching from differentiated to dedifferentiated melanoma cells induced by cytokines. Parameters are recapitulated in Table

••

  the first is performed on the ACT group and we focus on the dependence of the function g on the treatment dose D 70 . The optimal dose D the second optimization procedure is performed on the ACT+Re group and we focus on the dependence of the function g on the treatment dose D 70 , the restimulation dose D Re and the restimulation time t Re . The optimal treatment dose D opt * * 70 , retreatment dose D opt Re and retreatment time t opt Re are defined as follows: [D opt * * 70 , D opt Re , t opt Re ] = arg max (D 70 ,D Re ,t Re )∈D g(D 70 , D Re , t Re ), (4.5) where D = [D min , D max ] × [D Re min , D Re max ] × [t Re min , t Re max ]. We set the following initial values:t 0 = 0; n M 0 ∈ q α (n M 0 ,i ), α ∈ {5,50, 95}%, indicating different initial quantities of differentiated melanoma cells; D initial 70 = 0.02 the initial stimulation dose, D initial Re = 0.02 the initial restimulation dose and t initial Re = 160 the initial restimulation time. The treatment dose optimization interval [D min , D max ] is set to [0.005, 0.04] for both ACT and ACT+Re groups and the retreatment dose optimization interval [D Re min , D Re max ] is set to [0, 0.04] for ACT+Re group (such that D opt Re = 0 when D opt * * 70

Figure 4 .

 4 Figure 4. Individual fits (fitted tumor diameter (mm) along time (days)) for mice in CTRL group (a), ACT group (b), and ACT+Re group (c,d). Blue + marks represent no censored observations, red asterisk the censored observations, green line the fitted tumor size dynamics.

Figure 5 .

 5 Figure 5. Rescaled dynamics (population size/K) of differentiated cancer cells (top-left), dedifferentiated cancer cells (top-center), ratio of dedifferentiated over differentiated cells (topright), T cells (bottom-left), and cytokines (bottom-center). Scenarios: T cell survive (solid lines, population parameters of ACT group); T cell exhaustion (dotted lines, population parameters of ACT group and dT = to its quantile of order 99%). Vertical orange line: treatment start.

Figure 6 .

 6 Figure 6. T cell exhaustion probability (brown solid line) and global minimum of T cell population size (green dashed line) according to dT in ACT group (other parameters fixed to their estimated population mean). Black curve: T cell exhaustion probability as a function of dT . Red solid vertical lines: quantiles of order 5%, 50%, 95% of dT,i distribution; purple dotted vertical lines quantiles of order 1% and 99%.

)

  . Doses D70, DRe in number of T cells/K, restimulation time tRe in days.(a) ACT: q 95 (d M T,i ), q 95 (d T,i ), q 95 (l prod A,i ) (b) ACT: q 5 (d M T,i ), q 95 (d T,i ), q 50 (l prod A,i ) ACT+Re: D opt * * 70 = 0.04, q 95 (d M T,i ), q 95 (d T,i ), q 95 (l prod A,i ), q 95 (n M 0 ,i )

Figure 16 (

 16 Appendix) representing tumor size dynamics for different values of n M 0 supports this observation: the larger n M 0 , the larger tumor size at t = 70, the faster the tumor growth (even for D 70 = D opt * 70

Figure 8 .

 8 Figure 8. Individual fits. Blue + marks: no censored observations, red asterisk: censored observations, green fitted tumor size dynamics. From left to right, top to bottom: mice 1 to 5 in CTRL group, mice 6 to 12 in ACT group, and mice from 13 in ACT+Re group.

Figure 9 .

 9 Figure 9. Inidividual predicted values versus observed values. Blue points: no censored observations, red asterisk: censored observations.

Figure 10 .

 10 Figure 10. Visual Predictive Checks (pcVPC) for diagnosing the model. y-axis: tumor size (mm), x-axis: time (day). The blue and pink areas delineate respectively the confidence intervals for the 10th, 90th and 50th percentiles of the predictions obtained by simulations. The black curves represent the 10th, 50th and 90th theoretical percentiles. The green curves represent the empirical percentiles corresponding to the observed data and should remain within the corresponding prediction intervals to validate the model. The red areas indicate where the green empirical percentiles exit outside the corresponding prediction intervals (remember that the tumor size of the mice does not exceed 10 mm in the experiment). The blue and red points represent respectively uncensored and censored observations.

  Figure 12. T cell minimum according to treatment dose D70 in ACT group for dMT =

Figure 13

 13 Figure 13. T cell minimum according to treatment dose D70 in ACT group for dMT =

Figure 14 .

 14 Figure 14. T cell minimum according to treatment parameters in ACTRE group for dMT =

Figure 15 .

 15 Figure 15. T cell minimum according to treatment parameters in ACTRE group for dMT =

Figure 16 .

 16 Figure 16. Tumor size along time using quantiles q α (φi), α ∈ {5, 50, 95}%, φ ∈

  T trajectories of process T (t) and estimating the intermediate probability p k by considering the N S k trajectories which reach the intermediate threshold S k before time t F : pk =

	N S k

m. The S k (k < m) are the intermediate thresholds and S m = S the exhaustion threshold.

At each iteration k, the IS algorithm consists in simulating N

Table 1 .

 1 Estimated

		µ estimation ω estimation Quantile	Quantile
		(r.s.e (%))	(r.s.e (%))	5%	95 %
	r M	0.09 (4)	0.10 (34)
	r D	0.05 (10)	0.35 (20)
	s M D	< 0.01 (118) 0 (-)
	s DM -s M D < 0.01 (-)	0 (-)
	n M 0	0.08 (28)	0.58 (29)	0.03	0.21
	b T	< 0.01 (16)	0 (-)
	d T	0.02 (34)	1.07 (27)	< 0.01	0.09
	d M T	1.33 (55)	1.03 (48)	0.24	7.24
	s A	77 (37)	0 (-)
	d A	0.03 (-)	0.09 (181)
	l prod A	0.19 (93)	3.23 (25)	< 0.01	39.37
	σ	0.44 (6)
		Initial values		Optimisation range
	D 70	0.02		0.005	0.04
	D Re	0.02		0	0.04
	t Re	160		130	190

parameters for the final model: population mean µ and inter-mice standard deviation ω. r.s.e = standard error estimation

Table 2 .

 2 70 0.04 8.29 × 10 -7 6.77 × 10 -7 0.019 1.92 × 10 -5 4.45 × 10 -6 D opt Optimization results with quantiles q α (φi), α ∈ {5, 50, 95}%, and φ ∈

	q 5 , q 95 , q 95 , q 50 q 5 , q 95 , q 95 , q 95 q 50 , q 95 , q 95 , q 5 q 50 , q 95 , q 95 , q 50 q 50 , q 95 , q 95 , q 95 q 95 , q 95 , q 95 , q 5 q 95 , q 95 , q 95 , q 50 q 95 , q 95 , q 95 , q 95	Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re D opt 70 D opt Re t opt Re	--0.04 1.36 × 10 -6 1.10 × 10 -6 0.033 0.04 142 -0.04 -142 0.04 2.33 × 10 -6 1.86 × 10 -6 0.04 -0.012 -142 0.04 7.99 × 10 -7 6.61 × 10 -7 0.019 -0.04 -142 0.04 1.33 × 10 -6 1.08 × 10 -6 0.033 -0.04 -142 0.04 2.29 × 10 -6 1.83 × 10 -6 0.04 -0.019 -142 0.04 6.56 × 10 -7 5.78 × 10 -7 0.019 -0.04 -142 0.04 1.16 × 10 -6 9.84 × 10 -7 0.033 -0.04 -142 0.04 2.09 × 10 -6 1.72 × 10 -6 0.04 -0.019 -142	3.49 × 10 -5 4.73 × 10 -6 4.66 × 10 -5 5.22 × 10 -6 1.87 × 10 -5 4.44 × 10 -6 3.40 × 10 -5 4.72 × 10 -6 4.66 × 10 -5 5.20 × 10 -6 1.64 × 10 -5 4.38 × 10 -6 2.99 × 10 -5 4.65 × 10 -6 4.63 × 10 -5 5.13 × 10 -6
	{dMT , dT , l prod A column g initial , nM 0 } in ACT and ACT+Re groups. Column g opt ACT corresponds to g(D 70 ), opt * ACT to g(d initial ), column g opt ACT+Re to g(D opt * * 70 , D opt Re , t opt Re ) and g initial ACT+Re to
	g(d initial , d initial , t initial Re

i = ( i1 , . . . , in i ) ∼ N (0, σ 2 I n i ) represents the residual error η i = (η i r M , η i r D , η i s M D , η i s DM , η i n M 0

This work has been supported by the Malian Government and the LabEx PERSYVAL-Lab (ANR-11-61 LABX-0025-01). The authors would like to warmly thank Meri Rogava, Thomas Tüting and Michael Hölzel for fruitful discussions and for providing important medical informations, Jennifer Landsberg for providing the database used in this paper, and the Data Institute Univ.

A. Appendix

The code used to simulate the stochastic dynamics has been written in C++ and ran under the software R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] through the package Rcpp [START_REF] Eddelbuettel | Rcpp: Seamless r and c++ integration[END_REF].

Parameter (1/days)


1. Solar cell and module technology 2. Site selection: Albedo, Land cost, latitude… 3. Height, row distance, tracker, tilt, 4. Meteorological conditions A bifacial, single axis tracking (1T) pv plant can cost approximately 15% more than a comparable monofacial non tracking PV plant.

YIELD: a bifacial 1T plant can be well above 20%, up to 35 to 39% in some cases if conditions can be optimized. • High clearance of modules above the ground (up to several meters to give way to agricultural machinery)

BIFACIALITY IS HERE AND IT IS THE MAINSTREAM FOR

• Space between modules rows (to limit the shadow projected on the ground)

• Use of glass-glass modules for higher resistance (mist, agricultural chemicals, …)

 Bifacial modules are appropriate and may operate in good conditions to provide an energy gain.

System modelling is required for:

• Production estimations

  impact of potential-induced degradation (PID) -enabling longer lifespan. 2. No corrosion, no Aluminum 3. Easy to recycle 4. Less soiling and snow issues 5. increases the aesthetic appeal Drawbacks: 1. Specific structures 2. Cost higher / 2 glasses, high weight 3. edges of the panels can be little weak June 2022 / AgriVoltaics2022 conference Tandem technologies > 30% efficiency Maximizing solar spectrum use or varying which spectrum part is used for energy Silicon-based PV devices roadmap Efficiency, Cost, Reliability Sustainability, Eco-design, Recyclable Integrated, adapted • New paradigm for PV modules • Not more a standard for everything • Module has to be adapted to the application • The total cost of the system is calculated not only with LCOE basis (other KPI needed)• Reliability is a must • The system and O&M costs will depend

• System optimization • Operation monitoring PART 2 : MODELLING METHODS OF THE PV PRODUCTION

15

June 2022 / AgriVoltaics2022 conference Three kinds of models are needed:

• Thermal model • Electrical model • Optical model

	1.Large plants usually
	do not have rectangular
	shapes and show
	singularities (paths,
	missing arrays …)			4. Rebuilding the plant
				with the different types
		Real plant	Simulated plant	of zones
	2. Simulation of only a
	small elementary pattern with all edge effect possibilities:	Calculation of the irradiance on the module backside is a key issue: Simulated pattern	Conversion Table	3. Counting the different types of zones
	9 zones	• Mainly diffusion irradiance
		• Impact of neighbouring
		obstacles
		• Projected shadow on the
		ground
		• Irradiance mismatch
			June 2022 / AgriVoltaics2022 conference

June 2022 / AgriVoltaics2022 conference

June 2022 / AgriVoltaics2022 conference Solar generation capacity up to 25% depending on:


Morayo Adedjouma
email: morayo.adedjouma@cea.fr

Christophe Alix
email: christophe.alix@thalesgroup.com

Loic Cantat
email: loic.cantat@irt-systemx.fr

Eric Jenn

Juliette Mattioli
email: juliette.mattioli@thalesgroup.com

Boris Robert § †
email: boris.robert@irt-saintexupery.com

Fabien Tschirhart

Jean-Luc Voirin
email: jean-luc.voirin@fr.thalesgroup.com

Engineering Dependable AI Systems

Keywords: Artificial intelligence, Systems engineering and theory, System verification, Mission critical systems

If AI algorithms are now pervasive in our daily life, they essentially deliver non-critical services, i.e., services which failures remain socially and economically acceptable. In order to introduce those algorithms in critical systems, new engineering practices must be defined to give a justified trust in the capability of the system to deliver the intended services. In this paper, we give an overview of the approach that we have put in place to reach this goal in the framework of the French Confiance.ai program. Based on the needs of the industrial partners of the program, we propose a model-based analysis framework capturing the two dimensions of the problem: the one related to the development and operation of the system and the one related to the trust in the system.

I. INTRODUCTION

A. A long journey for deploying AI in an industrial solution

Artificial Intelligence (AI) faces challenges on different levels that need to be addressed to speed up the adoption and deployment of AI in industry. Towards that goal, we need a SoSE workbench and methods that support trustworthy AI at both the component and system levels, whether the specifications come from regulation, societal concerns, safety or security, etc.

With an allocated budget of C45M for the period 2021-2024, The French program "Confiance.ai" (www.confiance.ai) is a collective endeavour of 9 major French industrial (Air Liquide, Airbus, Atos, Naval Group, Renault, Safran, Sopra Steria, Thales and Valeo) and academic partners (CEA, IRT Saint Exupéry, IRT SystemX, INRIA). Its target is the creation of an engineering workbench enhanced by methods and tools that will enable the integration of AI into products or services considered as critical, i.e. where accidents, failures or errors could have serious consequences for people and/or valuable assets. The program is organized in 7 projects: EC1, which builds an integrated environment for the development and deployment of trustworthy AI (the "workbench" mentioned previously); EC2, focused on engineering processes and methods; EC3, focused on the characterization and qualification of AI components; EC4, focused on the design of AI components; EC5, focused on data, information and knowledge engineering; EC6, focused on on Integration, Verification, Validation and Qualification (IVVQ) activities; EC7, focused on embeddability issues. The results presented below have been achieved essentially in projects EC1, EC2 and EC6 in close collaboration with the other projects of the program.

Based on the observation that, in some areas, too few AI Proof of Concept (PoC) are reaching production level, a survey was sent to all industrial partners of Confiance.ai to identify their current practices and needs concerning trustworthy AI.

B. Capturing the industrial practices and needs

Each industrial partner received questionnaires dedicated to their respondent's business profile. Each question corresponded to a specific engineering phase as well as to an object of interest (algorithmic engineering, data engineering...) for a total of 140 different questions, most of them being open-ended. 52 questionnaires were returned, from which 3583 responses to individual questions were extracted. The following graph highlights the global distribution of answers according to the respondent's business profile. Organized in eighteen items, we started by trust and confidence issues, and general questions about the current and future use of AI solutions. Then, a number of objects of interest were addressed, namely the ODD (Operational Design Domain) specification, data, machine learning solutions, knowledge-based solutions, and some transverse aspects: human-computer interaction, safety and software interoperability. All engineering phases were examined, starting by lifecycle management, followed by seven phases: requirements, design, algorithmic engineering, assessment, deployment, supervision/monitoring and certification. The last item contained a few additional issues not included in the previous ones.

• Trust and confidence: Dealing with general considerations about trust/confidence, the questionnaire emphasizes the need to share a common view of the terminology used, of the available technologies and methods. However, most bibliographic references mentioned are norms and standards. The main expected services are linked to certification and validation on representative use cases, with supporting methods and tools. The collaboration between AI systems and human requests AI systems to be reliable, reactive, to reduce the worker's cognitive load, to be reproducible, transparent, in short to behave as a trustworthy co-worker. The major risk identified is that an AI system reproduces the same errors over and over again, leading to the famous saying "trust takes the stairs up and the elevator down". • Safety: The main safety concerns raised by industrial partners are the potential errors of an AI system and their consequences. Guaranteeing precision and generalization, characterizing the uncertainty of answers are obvious needs. Several standards have been designated, depending on the industrial domain and reveal that industrial partners almost systematically work in standard controlled environment.

• Lifecycle management: Lifecycle management is an activity that seems to be little addressed and, sometimes, relegated to a single tool (MLFlow and Git are often used for the model) and sometimes ignored. • Requirements engineering: Requirements are expressed in terms of performance, robustness, explainability, interpretability. When dealing with embedded systems, the usual constraints are mentioned: size, power consumption, real time performance. The methods and tools used for requirements engineering of AI components are the ones currently used by industrial partners for traditional components.

• Design: Formal methods, processes and methodologies and the supporting tools for AI components and AI-based systems are not often used in current practice. However, when companies integrate components from third parties (i.e. suppliers) they tend to ask for guarantees of some sort in the contracts. • Deployment: Deployment for embedded systems is the main concern since deployment on large-scale infrastructures or cloud does not seem to be a subject of concern. • Supervision/Monitoring: Generally speaking, the process of integrating the monitoring of AI-based systems does not seem mature. Various open-source tools are frequently mentioned (MLFlow, Kubernetes, Tensorboard, Prometheus, Grafana) as well as tools delivered with AI platforms (Azure, AWS Sagemaker), sometimes inhouse tools monitoring specific system indicators are mentioned. • Certification: It seems clear to everyone that certification of AI component will have a major impact on the current practices, for methodologies as well as for tools used in the development chain, essentially because there is little or no use of certification tools and methods for AI systems at the moment. Engineers are keen to apply such changes if the added-value is demonstrated. However, they ask to have some freedom of choice, not to be constrained by the requirements of certification tools and methods. All these issues highlight today's roadblocks that need to be addressed to bridge the gap between Proof of Concepts and actual deployment of dependable (system of) systems involving AI algorithms. We consider that those issues can only be tackled by a deep revamping of current engineering practices. The operational exploitation of AI is relatively recent. It was determined by the spectacular improvements of algorithms and the hardware components executing those algorithms. Initially exploited for non-critical tasks showing no or a very low level of risk, building an AI system was essentially a matter of combining ad-hoc engineering practices with the objective of providing "usable" results in the most cost-effective way. When it comes to industrial critical systems, several additional constraints must be considered. First processes must be rationalized, justified, made reproducible, optimized, etc. Second, processes must ensure that the overarching properties of the system under design are actually satisfied with the appropriate level of confidence: (i) the defined intended behavior of the system is correct and complete with respect to the desired behavior, (ii) the implementation of the system is correct with respect to its defined intended behavior, under foreseeable operating (iii) any part of the implementation that is not required by the defined intended behavior has no unacceptable safety impact [START_REF] Holloway | Understanding the Overarching Properties[END_REF].

In the context of Confiance.ai program, we propose to reconcile these two approaches, namely learning from ad-hoc engineering practices on the different use cases on the one hand, and structuring reproducible processes for achieving appropriate level of confidence on the other hand. To achieve this reconciliation, we need an analysis framework able to capture and organize engineering contexts, constraints, activities, data, lifecycle, etc. concurrently under different viewpoints, in order to build a global and comprehensive model. Each viewpoint enriches others in an iterative and incremental, multi-viewpoint analysis.

B. Capturing and modeling the engineering processes of dependable AI-based systems

To ensure that AI-based systems will possess the necessary "Trust Properties", specific System Development activities (elaborated by project EC2) and IVVQ activities (elaborated by project EC6) are required. These activities will use a chain of elementary methods and tools, specified or recommended by projects EC3, EC4, EC5 and EC7.

The roles and the relations between the main artefacts produced by these activities are depicted in Fig. 3 The approach consists in:

• Defining analysis viewpoints, formalized in a modeling tool by a meta-model containing the definition of involved concepts and semantic relationships between these concepts. As of today, the analysis framework implementation is derived from the Capella toolbox 1 . Capella has been chosen to develop the first version of the tool so as to leverage our experience on this technology and its underlying system development methodology (ARCADIA). However, our approach is independent from the modeling technology, and other tools may well be used. • Consolidating the methodological outputs of the Confiance.ai projects by analyzing their various aspects: engineering context, constraints, activities, data, lifecycle, etc. • Formalizing the analyzed methods in a modeling tool, according to the meta-models of the considered viewpoints. The modeling will help ensuring that all methods are compatible with each other, and that the Confiance.ai program will produce a consistent end-to-end process allowing the design of dependable AI-based systems. The overall structure of the analysis framework is depicted on Fig. 2.

The viewpoints for System Development (left part of the figure) are the following:

• 2 generic viewpoints:

-Engineering activities for trustable AI (layer 1): Define the tasks to perform so as to specify, design, produce, deploy and operate an appropriate and trustable solution to a well understood need, involving AI techniques. An example of engineering activities viewpoint limited to ML algorithm engineering is shown Fig. 4. -AI-related data life cycle (layer 2): Identify major data required/produced by AI engineering, when they are produced/used, and how they evolve with time

• 2 viewpoints dedicated to risk on trust (ie. risk that the trust on the capability of the system to deliver the expected service is reduced or lost):

-Risk on trust due to engineering (layer 3): identify major sources of bias or errors brought by other engineering activities to inputs and outputs of AI engineering and data -Risk on trust due to system during operation (layer 4): identify major sources of bias or corruption brought by other system components interacting with AI components in Operation

• 2 viewpoints dedicated to trust development and support: -Trust through system behavior (layer 5): define major system capabilities needed to ensure Trust in Operation -Trust through Human/System/AI Collaboration (layer 6): define expectations of humans stakeholders, their role and workshare with System AI, in delivering the expected services in a trustable manner. 4 transverse system viewpoints are also identified (right part of Fig. 2):

• Integration of AI functions: characterize and address specific concerns related to integrating one or more AI functions together in system target context; deliver guidance on how to manage each concern, • Trust on Performance, Safety, Security: define main needs, contributions and obstacles regarding Trust applied to AI decision performance, safety, and security of the global solution including AI, • Engineering Lifecycle Management: define processes to revisit engineering choices and decisions according to evolution of context, environment and needs, • Engineering Environment Capabilities: define the tooling support required to make trustable AI systems engineering feasible, scalable, efficient and secure. Consistency of the content in all those viewpoints shall be checked.

The next section focuses on VV aspects in relation with these engineering viewpoints.

III. COUPLING THE ENGINEERING AND V&V CONCERNS

Verification and validation activities are essential contributors to trust. Indeed, those means are essentially aimed at providing evidences that the system will realize the intended function. Towards that goal, we propose to establish a clear, traceable, auditable, and as formal as possible relationship between the engineering items produced by an AI system engineering activity, the properties that those items must satisfy, and the activities providing evidences that those properties are actually satisfied. In the Confiance.ai program, this relationship is captured by means of Assurance Cases [START_REF] Rushby | Assurance and Assurance Cases[END_REF].

An Assurance Case provides a structured argument to justify certain properties (sometimes called "claims") about the system under design, based on evidences concerning both the system and the environment in which it operates. The objective is to demonstrate as rigorously as possible that if some evidences are provided then some claim is justified. This argument cannot be as rigorous as a mathematical demonstration2 , simply because it does not refer to mathematical entities and mathematical properties. Nevertheless, the objective is to make it as close as possible to what a mathematical demonstration would be. In particular, terms and properties must be defined as precisely as possible, hypothesis and assumptions must be clearly stated, etc.

A claim concerns the satisfaction of some property by the system (e.g. system-level properties such as safety, security, or item-level properties such as "completeness", "consistency", etc.). Building an argument consists to decompose the initial claim into sub-claims deemed easier to justify in a divideand-conquer approach, down to the point where claims can be directly justified by showing evidences. During the construction of the argument, claims and properties may concern the whole system or some engineering items produced and used to engineer the system. As for the design of the system itself, building the Assurance Case of the system is not an ideal and strict top-down process that would start from some top-level property (e.g., "the system performs the intended function") and would be progressively decomposed into more and more primitive properties applicable to more and more primitive items. It is rather a combination of top-down and bottom-up approaches.

Assurance Cases are not a new practice. They actually derive from the now-well established practice of safety cases that are already required in some industrial domains (e.g. in the automotive domain with the ISO 26262 [START_REF]Road vehicles -Functional safety[END_REF]). For the interested reader, the genesis of Assurance Cases is very well described in [START_REF] Rushby | Understanding and evaluating assurance cases[END_REF].

As stated earlier, the introduction of AI in industrial practices strongly changes well-established practices, in particular those related to certification [START_REF] Mamalet | Machine Learning in Certified Systems[END_REF], and many initiatives are currently working on updating or defining new practices addressing the specificity of AI (e.g., EUROCAE WG114, ISO/TC204 It is worth noting that, in the context of Confiance.ai, Assurance Cases are strictly considered as a means to formalize the argumentation and build trust. It is a reference model from which other, specific representations, can be extracted. This includes, in particular, representations that will eventually be required by certification authorities. Our Assurance Cases are

• Not a verification and validation plan, but it may be used to build it. • Not a Certification Standard, but it may be used to build one and, at least, it may be used to determine the activities to be carried out to comply with the existing ones.

Traceability between the objectives / recommendations of standards and our Assurance Case is not yet addressed. The methods identified in the Assurance Case to provide evidences are normally captured by verification and validation activities of the engineering process. Therefore, the overall process is the following:

• The design / development / deployment / etc. workflow is established and the artefacts involved in this -workflow are described • Claims concerning the properties that those artefacts must possess are expressed • Assurance cases are built to show how those properties will be assessed • At the "bottom" of the argumentation (the leaves of the tree), one find evidences • Evidences are brought thanks to some dedicated verification and validation activity • Those activities are inserted in the workflow to give the complete picture. The navigation between the Assurance Case model and the engineering process model is illustrated on Figure 5. The Assurance Cases developed in Confiance.ai are fairly generic for they need to be applicable in different contexts (embedded system, production lines,. . . ), industrial domains (aeronautics, space, automotive, etc.), and for different types of applications involving different types of sensors and algorithms, with different level of criticality.

The Assurance Cases capture this diversity using "strategies". For instance, the same claim about a given property P may be decomposed (or justified) in different manners, or strategy, according to the level of criticality of the system. This way, a claim about temporal determinism may be achieved empirically by means of measurements performed on the actual system in one strategy, and may be performed using complex static analysis techniques based on formal methods (e.g., abstract interpretation) in another strategy. In this example, the confidence on the engineering item will (normally) depends on the strategy.

At the end of the day, the objective is to use the model to help making the optimal choice considering (i) the additional confidence brought by the method on the capability of the system to perform its intended function, and (ii) the cost of To reach this objective, we are currently developing a tool that allows navigating between the engineering workflow and the Assurance Case. A screenshot of the tool interface is given on Figure 6. It shows the workflow (on the top left), the Assurance Case (on the top right), and the applied strategy (on the bottom right). Thanks to this tool, the user will eventually be able to (i) build a V&V strategy considering the risk level associated with errors affecting engineering item, cost and trust indicators associated with the production of evidences, and to (ii) display the resulting engineering workflow including V&V activities.

IV. CONCLUSION: FROM THE MODEL TO THE WORKBENCH

Adoption of AI in our industry raises many technical and non technical challenges and a huge effort is currently deployed to address these challenges and facilitate the early industrial adoption of AI in a cost-effective and safe way.

In the context of the Confiance.ai program, the most critical of those challenges have been identified, covering two main aspects: dependability to provide the capability to give a justified trust on the capability of the system to deliver the expected service, and industrial efficiency in order to ensure that dependability will be achieved in a cost-effective way.

The Confiance.ai program addresses most of the dimensions of the problem, from the provision new ML algorithms and techniques addressing the various dimension of trust, including explicability, fairness, temporal determinism, etc. But the engineering practices themselves must also be updated to account for the specificities of AI. Towards that goal, we propose a "Trustable AI Engineering Definition Framework" to build system development and V&V workflow integrating explicitly the various dimension of trust. The framework relies on a model-based approach involving a series of 10 viewpoints capturing the various aspects of system development including those related to data engineering, risk analysis, etc.

As of today, a meta-model capturing and organizing the concepts supporting the different viewpoints has been developed, and its implementation in the Capella environment is on-going.

In the next phase, our objective is (i) to populate the engineering model using the methods and tools developed or recommended by the different projects (EC3,4,5,7), (ii) to build the complete V&V argumentation with clear links between the workflow activities, the engineering items, evidences, and activities providing evidences, and (iii) bring together related methodology elements and adequate tooling to support collaboration of all engineering disciplines for trustworthy AI-based products over their life cycle. The models obtained in this first phase will be used to evaluate, validate and possibly correct the approach.

Fig. 1 .

 1 Fig. 1. Global distribution of answers according to the respondent's business profile

Fig. 2 .

 2 Fig. 2. Global view of the analysis framework

  . The System Development Activities will produce Engineering Items. The IVVQ Activities will propose various Strategies to generate the Evidences showing that these Engineering Items possess the Trust Properties. The Trust Environment (produced by Confiance.ai project EC1) will orchestrate these System Development Activities and IVVQ Activities, and will store or reference the produced Engineering Items and Trust Evidences.

Fig. 3 .

 3 Fig. 3. Roles and relations between artefacts in Confiance.aiThe Trustable AI analysis framework designed by the Confiance.ai program will allow to elaborate the strategies for System Development Activities and for IVVQ Activities and will contribute to the specification of the Trust Environment.The approach consists in:

Fig. 4 .

 4 Fig. 4. Illustration of a 1st Capella model for ML Algo Engineering

Fig. 5 .

 5 Fig. 5. From the workflow to the Assurance Cases and backward

Fig. 6 .

 6 Fig. 6. Engineering and Assurance Case tooling implementing the method.To reach this objective, we are currently developing a tool that allows navigating between the engineering workflow and the Assurance Case. A screenshot of the tool interface is given on Figure6. It shows the workflow (on the top left), the Assurance Case (on the top right), and the applied strategy (on the bottom right). Thanks to this tool, the user will eventually be able to (i) build a V&V strategy considering the risk level associated with errors affecting engineering item, cost and trust indicators associated with the production of evidences, and to (ii) display the resulting engineering workflow including V&V activities.

See https://www.eclipse.org/capella

It is not strictly deductive.

This work has been supported by the French government under the "France 2030" program, as part of the SystemX Technological Research Institute