Chloé Pasin

François Dufour

Laura Villain

Huilong Zhang

Rodolphe Thiébaut
email: rodolphe.thiebaut@u-bordeaux.fr

Controlling IL-7 injections in HIV-infected patients

Keywords: Optimal control, immune therapy, dynamic programming

come

Introduction

The infection by the Human Immunodeficiency Virus (HIV) compromises the immune system functions, mainly because of the depletion of CD4 + T lymphocytes. Combined antiretroviral (cART) therapy has led to a spectacular improvement of patients' survival by controlling virus replication and consequently restoring the immune system functions. However, some patients fail at reconstituting their immune system and recovering normal CD4 + T cell levels, especially when they start antiretroviral treatment late [START_REF] Lange | Immune reconstitution with antiretroviral therapies in chronic hiv-1 infection[END_REF]). Immune therapy has been considered as a complement to cART to help immune restoration. In particular, interleukin-7 (IL-7), a cytokine produced by non-marrow-derived stromal and epithelial cells, is thought to improve thymic production [START_REF] Mackall | Il-7 increases both thymic-dependent and thymic-independent t-cell regeneration after bone marrow transplantation[END_REF]; [START_REF] Okamoto | Effects of exogenous interleukin-7 on human thymus function[END_REF]) and cell survival [START_REF] Tan | Il-7 is critical for homeostatic proliferation and survival of naive t cells[END_REF]; [START_REF] Vella | Cytokine-induced survival of activated t cells in vitro and in vivo[END_REF]; [START_REF] Leone | Increased cd4+ t cell levels during il-7 administration of antiretroviral therapy-treated simian immunodeficiency virus-positive macaques are not dependent on strong proliferative responses[END_REF]). The safety and beneficial effect of injections of exogenous IL-7 was first shown in phase I trials [START_REF] Sereti | Il-7 administration drives t cell-cycle entry and expansion in hiv-1 infection[END_REF]; [START_REF] Levy | Enhanced t cell recovery in hiv-1-infected adults through il-7 treatment[END_REF]) and observational studies [START_REF] Camargo | Responsiveness of t cells to interleukin-7 is associated with higher cd4+ t cell counts in hiv-1positive individuals with highly active antiretroviral therapy-induced viral load suppression[END_REF]). Then, phase I/II human clinical trials (INSPIRE 1, 2 and 3 studies) have evaluated the effect of repeated cycles of three IL-7 injections and showed that this therapy helped maintaining HIV infected patients with CD4 + T cells levels above 500 cells/µL [START_REF] Levy | Effects of recombinant human interleukin 7 on t-cell recovery and thymic output in hiv-infected patients receiving antiretroviral therapy: results of a phase i/iia randomized, placebocontrolled, multicenter study[END_REF]), a level associated with a nearly healthy clinical status [START_REF] Lewden | Hiv-infected adults with a cd4 cell count greater than 500 cells/mm3 on long-term combination antiretroviral therapy reach same mortality rates as the general population[END_REF]).

The dynamics of CD4 + T lymphocytes following IL-7 injections can be fitted by mechanistic models based on ordinary differential equations (ODEs). These models contain compartments corresponding to different populations of CD4 + T lymphocytes and biological parameters characterizing these populations. Hence it was possible to quantify the effect of repeated cycles of IL-7 on CD4 + T lymphocytes on specific parameters. Previous work using data from clinical trials (INSPIRE studies) has shown that IL-7 enhances both proliferation and survival of CD4 + T lymphocytes [START_REF] Thiebaut | Quantifying and predicting the effect of exogenous interleukin-7 on cd4+ t cells in hiv-1 infection[END_REF]). Moreover, a differential effect of the injections within a given cycle has been found, the third injection of a cycle appearing to have a weaker effect on proliferation than the first ones [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF]).

In addition to providing insight into the most important mechanism of the effect of exogenous IL-7, the models have shown a very good predictive capacity [START_REF] Thiebaut | Quantifying and predicting the effect of exogenous interleukin-7 on cd4+ t cells in hiv-1 infection[END_REF]; [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF]). Hence, the next step was the determination of the best protocol of injections. A first approach, realized in [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF], consisted in simulating and comparing the regular protocol to three other protocols with different numbers of injections by cycle.

In all four protocols, CD4 counts were measured every three months, and a new cycle was administered when the CD4 numbers were below 550 cells/µL. Comparison was based on three criteria: number of injections received, mean CD4 count and time spent below 500 cells/µL over a four-year period. Results showed that cycles of two injections could be sufficient to maintain CD4 levels, while using less injections than in the clinical protocol. These results suggest the possibility to reduce the number of injections in clinical protocols. However, the three months delay between visits is independent of the patient and constrains the protocol. While some patients with "not too low" baseline CD4 levels could afford coming back later than three months after the last visit, some patients with "low" baseline CD4 levels would need more repeated cycles or more injections by cycle. Individualized protocols could help in achieving the maintenance of the patient's CD4 + T lymphocytes levels over a given threshold by using different patient-dependent timing of injections and doses. The possibility of conducting the lightest intervention for every patient could be very important for the development of IL-7 in HIV infected patients especially for further large clinical trials.

Optimization of schedule and doses is an up-to-date question when working on protocol of injections. In their review on mathematical modeling for immunology, Eftimie et al (2016) emphasize the need for complex optimal control approaches coupled with immunology experiments, in order to improve clinical interventions. Basically, there are two kinds of techniques that can be used to solve optimal control problems: methods involving Pontryagin's maximum principle and dynamic programming approaches. Pontryagin's maximum principle has been applied to a number of biological problems of the form dx(t) dt = f (x(t), u(t)), where the solution to the ordinary differential equation depends on the dynamics of the control function u(t). For example, it was applied to the determination of the optimal schedule of dendritic cells vaccine injection in cancer immunotherapy by [START_REF] Castiglione | Optimal control in a model of dendritic cell transfection cancer immunotherapy[END_REF]; [START_REF] Cappuccio | Determination of the optimal therapeutic protocols in cancer immunotherapy[END_REF]; [START_REF] Castiglione | Cancer immunotherapy, mathematical modeling and optimal control[END_REF] and [START_REF] Pappalardo | Vaccine protocols optimization: in silico experiences[END_REF]. However, in our case, the model is a piecewise deterministic Markov model (PDMP), where dynamics of IL-7 are unknown and not modeled. Addressing the objective of spending the least time possible under the threshold of 500 cells/µL by using repeated injections of IL-7 corresponds in a more formal way to determining actions (injection or not and choice of dose) at given time points over a horizon of time: this can be treated as a problem of impulse control in the optimal control theory. To the best of our knowledge, there is no maximum principle solving this kind of problem. We will focus on a dynamic programming method, as developed in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF]. In a formal mathematical framework, we addressed the question of optimizing the schedule of IL-7 injections for a given patient by a two-steps method: determining an adapted mathematical model for the process, and developing a numerical method to determine an optimal strategy of IL-7 injections for a given patient.

As described in [START_REF] Davis | Piecewise-deterministic markov processes: A general class of non-diffusion stochastic models[END_REF], most of the continuous-time stochastic problems of applied probability (including those modeling biological processes) consist of some combination of diffusion, deterministic motion and/or random jumps. Ordinary differential equations can be included in the class of deterministic motion with random jumps. In our particular framework of modeling cell dynamics after IL-7 injections, jumps correspond to the change of some parameters value. This can be easily and naturally modeled by the largely studied class of Piecewise Deterministic Markov Processes (PDMPs). A non-controlled version of this model can be described by iteration as follows: from a point in the state space, the process follows a deterministic trajectory determined by the flow, until a jump occurs. This jump happens either spontaneously in a random manner, or when the flow hits the boundary of the state space. After the jump, the system restarts from a new point determined by the transition measure of the process. We will show in this article how to model the dynamics of the CD4 + T cells in HIV-infected patients following IL-7 injections using a PDMP.

According to the problem studied in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF], impulse control consists in possible actions only when the process reaches its boundary. This will constitute our framework: the decision maker has the possibility to inject IL-7 when the number of CD4 + T lymphocytes reaches a given level or when a certain amount of time has passed since the last injection. Each action has a cost, and a strategy is defined as the set of all realized actions over a given horizon. The impulse problem consists in determining a strategy of injections minimizing the optimality criterion induced by the cost function. In our case, the cost function depends on the number of injection realized and the time spent with the CD4 + T lymphocytes levels under the threshold of 500 cells/µL, as both quantities should be minimized.

As emphasized by the authors of [START_REF] Dufour | Numerical methods for simulation and optimization of piecewise deterministic markov processes[END_REF], the development of computational methods for the control of PDMPs has been limited, and at the moment there is no general method allowing the numerical resolution of optimal control on PDMPs (and in particular impulse control). This constitutes a real challenge. We propose in this work a numerical method based on the results developed in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF]. In this paper, the authors studied the existence of a solution of the Bellman-Hamilton-Jacobi equation by showing that the value function is the limit of a sequence of functions given by iteration of an integro-differential operator. This construction leads to a natural method for the computation of the optimal cost and the determination of an optimal strategy of injections. In particular, we have developed an algorithm for the iteration of the operator and applied our numerical tool to the case of the biological model. This provides a proof of concept as it succeeded in determining an optimal strategy for a number of pseudo-patients simulated using previous estimations. The paper is organized as follows: section 2 presents the mathematical modeling of the process, including data and design of INSPIRE studies, as well as mechanistic model and finally the associated PDMP. Section 3 focuses on the optimal control problem, by reminding the main theoretical results from [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF] and adapting them to the IL-7 study. Section 4 presents some numerical aspects of the dynamic programming work, necessary to determine the optimal cost function and strategy for a given patient. Results are presented in section 5 and discussion is done in section 6.

Mathematical modeling

Material

Our work is based on three phase I/II multicenter studies assessing the effect of a purified glycosylated recombinant human Interleukin 7 (IL-7) treatment for immune restoration in HIV infected patients under treatment: INSPIRE [START_REF] Levy | Effects of recombinant human interleukin 7 on t-cell recovery and thymic output in hiv-infected patients receiving antiretroviral therapy: results of a phase i/iia randomized, placebocontrolled, multicenter study[END_REF]), INSPIRE 2 and INSPIRE 3 [START_REF] Thiébaut | Repeated cycles of recombinant human interleukin 7 in hiv-infected patients with low cd4 t-cell reconstitution on antiretroviral therapy: results of 2 phase ii multicenter studies[END_REF]). A total of 128 HIV-infected patients under antiretroviral therapy with CD4 + T cell count between 100-400 cells/µL and undetectable viral load for at least 6 months were included among the three studies from the time of the first injection. IL-7 was administered in cycles of weekly injections, with a "complete cycle" defined as three weekly injections. In INSPIRE, all 21 patients received complete cycles of IL-7 at different weight-dependent doses: 10, 20 and 30 µg/kg. In INSPIRE 2 and INSPIRE 3, 23 and 84 patients (respectively) received repeated (and sometimes incomplete) cycles of IL-7 at dose 20 µg/kg. Repeated visits and follow-up once every 3 months after the first cycle allowed to measure biomarkers levels in patients, in particular total CD4 + T cell counts and number of proliferating CD4 + T cells through Ki67 marker. At every visit, a new cycle of injections was administered if the patient's CD4 + T cell level was under 550/µL, in order to globally maintain the levels above 500 cells/µL. The total duration of the studies was 12, 24 and 21 months for INSPIRE, INSPIRE 2 and INSPIRE 3, respectively.

Mechanistic model

The dynamics of CD4 + T lymphocytes were largely described in [START_REF] Thiebaut | Quantifying and predicting the effect of exogenous interleukin-7 on cd4+ t cells in hiv-1 infection[END_REF] and [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF] by using several mechanistic models. We focus here on the following model, described in figure 1: it includes two populations of cells, non-proliferating (or resting, R) and proliferating (P). Resting cells are produced by thymic output at rate λ, become proliferating cells at rate π and die at rate µ R . Proliferating cells die at rate µ P and can also divide and produce two non-proliferating cells at rate ρ. The system of differential equations is written as follows :

     dR dt = λ -µ R R -πR + 2ρP dP dt = -µ P P -ρP + πR (1)
We assume the system is at equilibrium at t = 0, before the study begins and any injection is administered. IL-7 injections are realized through cycles containing up to three injections with seven days elapsed between each injection. Parameters estimation was performed using a population approach. Mixed-effect models including intercept, random and fixed effects, were used on log-transformed parameters, in order to both obtain an estimation across population and account for between-individuals variability. In the controlled framework, the decision-maker can decide to inject IL-7 to a patient at a given dose d, and this will affect the value of the proliferation rate π. Each injection denoted by n ∈ {1, 2, 3} of a cycle has a different effect on the value of π for patient i, defined as follows:

R P λ µ R µ P π ρ
πi = π0 + β (n) π d 0.25 i 1 {t∈[t i inj ,t i inj +τ i ]} (2)
with π = log(π);

β π = [β (1) π , β (2)
π , β

π ] is the vector of effect of each injection of a single cycle; d is the injected dose; t inj is the time (in days) at which IL-7 is injected and τ the length of effect of the injection (in number of days), considered equal to 7 in previous models [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF]. Estimation of parameters showed that effect of successive injections on the proliferation rate decreases within a cycle, and in particular the third injection seems to have a much weaker effect (as [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF]).

β (3) π < β (2) π < β (1) π )

Mathematical model: piecewise deterministic Markov process

As described in the introduction, ODEs-based mechanistic models can be included into the broader class of PDMPs. A PDMP is characterized by a state space in which it evolves, a flow, a jump intensity and a measure of transition. From a mathematical point of view, we note X the state space, an open subset of R d , d ∈ N, and ∂X its boundary. The flow associated with the process is φ(x, t) : R d × R → R d . The active boundary is defined as Ξ = {x ∈ ∂X : x = φ(y, t)} for some y ∈ X and t ∈ R * + . We will then denote X = X ∪ Ξ and for x ∈ X, we can define t * (x) = inf{t ∈ R + : φ(x, t) ∈ Ξ}. The controlled jump intensity η is a R + -valued measurable function and determines the law of the stochastic jumps. When the process, i.e. the trajectory of CD4 + T lymphocytes, reaches Ξ, the decision-maker can act by injecting IL-7 to the patient. The action varies according to the dose injected. This leads to a jump in some parameters value, and the process restarts from a new point defined by the transition measure Q(.|φ(x 0 , τ ), d), depending on the dose and the position of the state before the jump φ(x 0 , τ ).

In this section, we present the PDMP associated to the biological process described in section 2.2. Here, the PDMP is patient-dependent. As we focus on the control question (and not the estimation one), we suppose that parameters values of the studied patient are known. Previous work has shown that estimation of ODE's parameters based on population approaches can be reliable [START_REF] Lavielle | Estimation of population pharmacokinetic parameters of saquinavir in hiv patients with the monolix software[END_REF], [START_REF] Prague | Nimrod: A program for inference via a normal approximation of the posterior in models with random effects based on ordinary differential equations[END_REF]). Moreover, the model developed in our particular framework for the CD4 dynamics has shown good predictive abilities [START_REF] Jarne | Modeling cd4 + t cells dynamics in hiv-infected patients receiving repeated cycles of exogenous interleukin 7[END_REF]). Therefore we make the assumption that we determine the strategy for a patient who is already included in a clinical study and for which we had enough observations to estimate its parameters (by running a first cycle of injection for example). As developed in this part, the stochasticity is supposed to be induced by the biological model but not by the uncertainty on the parameters estimation. Sensitivity analysis of the method regarding the estimation uncertainty is provided in Appendix C.

The PDMP modeling the dynamics of CD4 + T lymphocytes of a given patient is defined using six variables: the state vector is denoted by x = (γ, n, σ, θ, p, r). γ determines the value of parameter π when combined with n, the number of injections realized in the ongoing cycle.

If d = [d 0 , d 1 , .., d m d ] is
the vector of all possibles doses (with d 0 = 0), then γ ∈ {1..m d + 1}. Injecting dose d k at the n-th injection of a cycle gives the following: 0.25 . The two variables σ and θ are time variables, discretized with steps of one day. In particular, σ corresponds to the number of days since the last injection and θ to the running time (θ = 1 at the first injection of the first cycle). Finally, variables p and r are values of compartments P and R solutions of system 1 with parameter π defined by γ and n and other parameters are supposed to have been previously estimated. We suppose the patient is followed until a horizon of time T h , then the state space is

γ(d k ) = k + 1 and π = π 0 +β (n) π d(γ)
X = X ∪ ∆ with X = {1..m d + 1} × {1..n inj } × {0..T h -7(n inj -1)} × {0..T h } ×{p min ..p max } × {r min ..r max }
and ∆ is an absorbing state representing the end of the study, at t = T h : ∆ = (0, 0, 0, T h , 0, 0). For x = (γ, n, σ, θ, p, r) ∈ X, the flow is defined as :

-φ(x, t) = (γ, n, σ + t, θ + t, p, r) if θ ≤ 1 -φ(x, t) = (γ, n, σ + t, θ + t, P (t, γ, n), R(t, γ, n)) if θ ∈ [1, T h -1]
, with P (t, γ, n) and R(t, γ, n) solutions of system 1 with initial conditions p and r and π determined with γ and n -φ(∆, t) = ∆ Moreover, even if the deterministic mechanistic model allowed good fits for the data, we make the hypothesis that the process undergoes some stochasticity: in particular, as the value of parameter π is modified by an injection of IL-7 during some days, we suppose that this modification can last randomly up to 7 days after the injection. Stochastic jumps can then occur with intensity η such that for x ∈ X, η(x) = η1 {γ>1} with η a given value and η(∆) = 0. It means that if we consider the modification of π value after an IL-7 injection through equation 2, τ follows there a random exponential law of parameter η. We define the constant K = η such that η(x) ≤ K for every x ∈ X. IL-7 injections aim at maintaining the CD4 + T cell level over 500 cells/µL. When this value is reached, we consider that the system has reached a deterministic boundary of the state space. A new injection of IL-7 injection is possible at that moment and gives the possibility to increase CD4 + T cell counts. To account for clinical constraints, we assume a minimum time σ min is observed between the beginning of two consecutive cycles, even if the number of CD4 falls below the threshold of 500. During cycles, the deterministic boundary corresponds to the seven days delay between injections. In a more formal way, the boundary can actually be reached in five different situations described in the following :

-for a technical reason due to the mathematical modeling which cannot account for an impulse action at t = 0, we define a first artificial boundary when the study begins, at θ = 1 : Ξ 1 = {x : θ = 1}. This allows a cycle of injections to begin at θ = 1. We suppose the studied patient is already included in the clinical study: it means that its biological parameters are known, and her/his CD4 + T cell count at t = 0 as well (either because she/he is at equilibrium, and the values are known from biological parameters, or because some measures have been realized at this time). -we also define a time corresponding to the end of the study and a boundary when the time reaches the horizon T h :Ξ 2 = {x : θ ≥ T h } -another boundary is reached when the patient is undergoing a cycle of injections and seven days have passed since the last injection : Ξ 3 = {x : n < n inj , σ = 7, θ < T h } -we also consider a boundary when at least one cycle was already achieved and the count of cells is equal to or below the threshold of 500 cells/µL. We also assume a minimum time σ min is observed between the beginning of two consecutive cycles : Ξ 4 = {x : p + r ≤ 500, n = n inj , σ ≥ σ min , θ < T h } -finally, an artificial boundary is created when π has not returned to its baseline value seven days after the last injection of a cycle :

Ξ 5 = {x : γ > 1, n = n inj , σ = 7, θ < T h }
We define the active boundary as

Ξ = Ξ 1 ∪ Ξ 2 ∪ Ξ 3 ∪ Ξ 4 ∪ Ξ 5 .
In this process, actions (IL-7 injections) can only be realized when the process hits the active boundary. We model the possibility of not doing an injection in a given cycle by using a fictive dose d 0 equal to zero. When beginning a new cycle of injections, the first injection needs to be positive though. The possible action made by the decision-maker depends on the boundary reached. Therefore for every x ∈ Ξ :

A(x) =    {d 1 , ..d m d } if x ∈ Ξ 1 ∪ Ξ 4 {0, d 1 , ..d m d } if x ∈ Ξ 3 ∅ if x ∈ Ξ 2 ∪ Ξ 5
We also define the transition measure (or Kernel): it determines the new point from which the process restarts after a jump. It depends on the injected dose only when the boundary of the process is reached. All possible situations are the following :

-when the flow hits Ξ 1 , the study begins with administration of a cycle of injections. γ takes the value corresponding to the chosen dose. (p, r) = (P c , R c ), known values from either equilibrium or biological measures made on the patient before the beginning of the study -when the flow hits Ξ 2 , the study is over and nothing happens from absorbing state ∆ -when the flow hits Ξ 3 , a new injection is administered to the patient. γ takes the value corresponding to the chosen dose γ(d), n increases by one, σ goes back to 0 -when the flow hits Ξ 4 , a new cycle of injections begins. γ takes the value corresponding to the chosen dose, n goes back to 1, σ goes back to 0 -when the flow hits Ξ 5 , there is no injection. γ goes back to 1 -in case of spontaneous jump, there is no injection and γ goes back to 1

In a formal way, the Kernel Q is written :

Q(dy|x, d) = δ (γ(d),1,0,1,Pc,Rc) (dy)1 {x∈Ξ1} + δ ∆ (dy)1 {x∈Ξ2} + δ (γ(d),n+1,0,θ,p,r) (dy)1 {x∈Ξ3} + δ (γ(d),1,0,θ,p,r) (dy)1 {x∈Ξ4} + δ (1,n,σ,θ,p,r) (dy)1 {x∈Ξ5} + δ (1,n,σ,θ,p,r) (dy)1 {x∈ X}
The impulse control problem consists in determining the optimal scheme of injections and their associated dose according to a given optimality criterion, based on the cost function C : in our case, this cost function depends on the number of injections realized and the time spent with the CD4 + T lymphocytes levels under the threshold of 500 cells/µL. Both quantities need to be minimized, in order to maintain the patient in good health by injecting the least possible. The cost can be divided in two parts. First, the gradual cost penalizes the trajectory of the process through the time spent under the threshold after the beginning of the first cycle. This time is considered in months, approximately, as it is computed as the number of days divided by 30. For x = (γ, n, σ, θ, p, r) ∈ X :

C g (x) = 1 30 1 {p+r<500} 1 {θ≥1}
Then, the cost associated to an impulsive action penalizes the fact of injecting IL-7 to the patient :

C i (x, d) = 1 {x∈Ξ1∪Ξ4} + 1 {d =0} 1 {x∈Ξ3}
After the horizon, the cost is null, as C i (∆) = C g (∆) = 0.

Optimal control

In this section, we will first remind the main theoretical results obtained in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF], then we will transpose these results to our particular context.

Main theoretical results

The objective of this section is to adapt some results obtained in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF] to our specific context. We follow closely their notation. The set of all realized injections over a given horizon constitutes the strategy of injections.

In a more formal way, a strategy u of the decision maker is a sequence u = {u n } n∈N of functions u n : X → A giving the action to realize at punctual time points t n ≥ 0 when the system is in state x ∈ X. The set of all admissible strategies is noted U. According to section 2.2 in Costa et al ( 2016) there exists a continuous-time stochastic process ξ defined on probabilistic space using characteristics φ,η and Q depending on the action given by u, such that ξ t , t ∈ R + corresponds to the state of the variables at time t. To each admissible strategy u ∈ U, we associate a discounted cost optimality criterion depending on the gradual cost on the trajectory of the process ξ, C g , and the cost related to an injection, C i , as defined in section 2.3:

V(u, x 0 ) = E u x0 ]0,+∞[ e -αs C g (s)ds +E u x0 ]0,+∞[ e -αs I {ξs-∈Ξ} A(ξs-) C i (ξ s-, a)u(da|s)µ(ds) (3)
with α > 0 the discount factor and where µ is a measure that counts the number of jumps in the process. The impulse control problem aims at finding a strategy u minimizing the discounted cost optimality criterion. Here we want to determine the patient-specific schedule of injections and their dose to optimize the patient's CD4 + T lymphocyte numbers by using a minimum number of injections. The theorem allowing to determine the optimal cost and providing an optimal strategy is adapted from theorem 5.5 in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF]. It is stated as followed:

Theorem 1 Suppose assumptions A, B and C from section 3.2 in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF] are verified. We define the sequence of functions {W q } q∈N for any x ∈ X as follows:

W q+1 (x) = BW q (x) for q ∈ N W 0 (x) = -K A 1 Aε 1 (x) -(K A + K B )1 A c ε 1 (x) (4)
with constants K A and K B defined as in section 5 of [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF], K+α)t RV (φ(y, t))dt + e -(K+α)t * (y) TV (φ(y, t * (y))) ( 5)

A ε1 = {x ∈ X : t * (x) > ε 1 } and BV (y) = [0,t * (y)[ e -(
with real-value functions RV and TV defined for any V respectively on X and Ξ :

RV (x) = C g (x) + qV (x) + ηV (x) TV (z) = inf d∈A(z) C i (z, d) + QV (z, d)
q being the signed kernel, which computes the difference between the states before and after the spontaneous jump. For x ∈ X, it is defined with :

q(dy|x) = η(x)[Q(dy|x) -δ x (dy)]
The sequence of functions {W q } q∈N converges to a function W defined on the state space and such that : i) W (x 0 ) = inf u∈U V(u, x 0 ), with V defined as in equation 3 ii) there is a measurable mapping ϕ : Ξ → A such that ϕ(z) ∈ A(z) for any z ∈ Ξ and satisfying

C i (z, ϕ(z)) + QW (z, ϕ(z)) = inf d∈A(z) C i (z, d) + QW (z, d) . (6)
This theorem allows to determine the optimal cost and an optimal injection strategy, consisting in choosing the optimal action ϕ(z) for every point z ∈ Ξ reached on the trajectory of the process. Indeed, the iteration of the sequence {W q } q∈N defined by equation 4 can be realized by numerically approximating the operator B defined in equation 5. This will give an approximation of the function W , and in particular of W (x 0 ), corresponding to the optimal cost. Moreover, to obtain an optimal strategy, the process is the following: we simulate a trajectory from x 0 , then when a boundary is reached, the chosen action corresponds to the one minimizing the criterion C i (z, d) + QW (z, d), as given by equation 6.

Application

The process describing the effect of IL-7 on CD4 + T lymphocytes dynamics is now well defined by its characteristics φ, η and Q, boundaries and possible actions in section 2.3. Moreover, both gradual cost on the trajectory and impulse cost were defined in that section. We will quickly describe in this part how to apply the results from theorem 1 for our specific problem, i.e. determining the function B needed for the computation of the optimal strategy. For a more detailed and formal computation, we refer the reader to appendix A. We need to compute: K+α)t RV (φ(y, t))dt + e -(K+α)t * (y) TV (φ(y, t * (y)))

BV (y) = [0,t * (y)[ e -(
We define

G(V, y) = [0,t * (y)[ e -(K+α)t RV (φ(y, t))dt (7)
and

H(V, y) = e -(K+α)t * (y) TV (φ(y, t * (y)) (8)
We define a time interval ∆t (in practice equal to one day) and for every y = (γ, n, σ, θ, p, r) ∈ X, we note

n * (y) = t * (y) ∆t
For every j ∈ {0..n * (y) -1}, we denote φ j (y, t) = φ(y, j∆t) and φ(y, t * (y)) = (γ, n, σ + t * (y), θ + t * (y), p * (y), r * (y)). The integral defined in equation 7 is computed by approximation using the classic trapezoidal rule using the j∆t nodes. Thus, G(V, y) can be approximated by a linear combination of {V (y j )} j∈{0..n * (y)-1} , with y j depending on φ j (y, t). Moreover, H(V, y) is proportional to V (y), with y depending on the boundary reached in φ(y, t * (y)).

Finally, for every point y ∈ X, if we note y = y n * (y) , BV (y) can be computed as a linear combination of {V (y j )} j∈{0..n * (y)} .

4 Numerical aspects of the dynamic programming method

From theorem 1 we know that we need to compute the sequence {W q } q∈N such that for y ∈ X, W 0 (y

) = -K A 1 Aε 1 (y) -(K A + K B )1 A c ε 1 (y)
and W q+1 (y) = BW q (y) for q ∈ N. The sequence converges to a function W defined on X that allows the determination of the optimal cost and the optimal protocol of injections achieving that cost. This computation is realized on a grid of the state space: at each iteration q, a new matrix is computed, each element on line v and column s corresponding to BW q (x vs ), with x vs element of the grid Γ of the state space. The implementation of our algorithm was realized in Matlab version R2016b (The MathWorks, Inc., Natick MA, USA, 1984). In this section, we give elements to understand how the method is implemented. The structure of the code is detailed in Appendix B.

Discretization of the state space

The grid Γ contains points of the form (γ, n, σ, θ, p, r). γ and n are discrete variables with γ ∈ {1...m d + 1}, n ∈ {1..n inj }. σ and θ are discretized with a time step of one day, with σ ∈ {0..σ max } and θ ∈ {0..T h }. Solutions p and r of the ODE system are continuous and are discretized in a regular grid, with p ∈ {p min ..p max } with regular step h p and r ∈ {r min ..r max } with regular step h r . We then obtain :

n p = p max -p min h p + 1 n r = r max -r min h r + 1
h r and h p are chosen such that both n p , n r ∈ N count the number of values of p and r on the grid, respectively.

Organization of the grid

We arrange all points of the grid Γ in a matrix M of size N sum × N pr , with N sum corresponding to the number of possible (γ, n, σ, θ) combinations and N pr = n p n r number of possible (p, r) combinations. Each element M (v, s) v∈{1..Nsum} s∈{1..Npr} corresponds to a given combination (γ, n, σ, θ, p, r) of Γ , through the following bijection :

χ : Γ → {1..N sum } × {1..N pr } x vs = (γ, n, σ, θ, p, r) → (v, s) = χ l (γ, n, σ, θ), χ c (p, r)
χ l is defined in the following way : v corresponds to a given value of (γ, n, σ, θ).

Possible combinations of (σ, θ) depend on the value of (γ, n): for example, during the first cycle, when n = 1, σ = 0 is associated with θ = 1, while when n = 2, σ = 0 is associated with θ = 8. We divide the lines of matrix M by defining then N γn = (m d + 1)n inj blocks, corresponding to the possible combinations of (γ, n). Each block is indexed by i

= f (γ, n) = γ +(m d +1)(n- 1) ∈ {1.
.N γn } and contains combinations of (σ, θ), indexed by j = g i (σ, θ) ∈ {1..N bi } within the i-th block. The total number of lines of matrix M is the sum of the number of lines in each block: N sum = Nγn i=1 N bi . We can define a vector l block = (1, 1 + N b1 , .., 1 + k i=1 N bi , ..,

Nγn-1 i=1

N bi ) of length N γn , that determines the index of the first line of each block. Finally :

v = χ l (γ, n, σ, θ) = l block (i) + j -1
with i = f (γ, n) and j = g i (σ, θ). χ c is defined in the following way :

s = χ c (p, r) = p -p min h p + 1 + n p r -r min h r
such that s = 1 when (p, r) = (p min , r min ) and s = n p n r when (p, r) = (p max , r max ).

Iteration of the algorithm

Each iteration of the algorithm computes then a matrix M q such that M q (v, s) = W q M (v, s) = W q (x vs )

For every x = (γ, n, σ, θ, p, r) ∈ Γ , W q (x) is a linear combination of some W q (x m ), m ∈ {1..M x }, as shown in equation 12 from appendix A. Values of W q (x m ) are given by M q (χ(x m )); they are linearly combined and implemented in M q+1 (χ(x)).

Convergence criterium

We assume that the sequence converges when W q+1 -W q ∞ < . In practice, we compute max v,s |M q+1 (v, s) -M q (v, s)| and we consider that the sequence converges with = 0.001. It usually occurs after 35 to 45 iterations.

Results

We applied the previously described method to the model detailed in section 2.3, with a choice of m d = 2 possibles doses: d = [0, 10, 20] (unit = µg/kg), cycles of 3 injections: n inj = 3 and a reduced horizon T h = 365 days. We also assumed a minimum time of σ min = 30 days between the end of a cycle and the beginning of a new one. For a given patient with fixed biological parameters, we can approximate the function W in a grid of the state space through convergence of the sequence {W q }: this determines the optimal cost over all strategies. Moreover, using equation 6 from theorem 1, we can simulate the strategy choosing the optimal action to realize when reaching the boundary of the state space and compute the cost of the obtained strategy. As some randomness is included in the model by the time of effect of an injection of parameter π, we simulate N = 5000 realizations of a protocol on a given patient with a Monte Carlo method and compute the expectation of its cost. From that, we check the numeric performance of our method by first comparing the cost of the optimal strategy to the computation of the optimal cost from the value function W . Moreover, we wish to compare the optimal strategy to other "naive" protocols. For each protocol, including the optimal one, we compute by Monte Carlo the mean cost, the standard deviation and the minimum cost achievable. This is usually reached when the patient responds well to all injections, i.e. the effect of the injection on parameter π lasts 7 days after every injection. In order to compare protocols based on clinical criteria, we also computed by Monte Carlo the mean number of CD4+ T cells count until horizon, the mean time spent under 500 cells/µL (in days) and the mean number of injections over all simulations. These comparisons were realized with 50 pseudo-patients. Parameters values were generated from the posterior law estimated on real data from INSPIRE trials in [START_REF] Thiébaut | Repeated cycles of recombinant human interleukin 7 in hiv-infected patients with low cd4 t-cell reconstitution on antiretroviral therapy: results of 2 phase ii multicenter studies[END_REF]. Patients are divided in three categories according to their initial levels of CD4 + T cells: "very low" baseline (100 -200 cells/µL), "low" baseline (200 -300 cells/µL) and "not too low" baseline (300 -400 cells/µL). Table 1 sums up the characteristics of the pseudo-patients population.

We first compare the value of the optimal function obtained from the numerical computation of W with the cost of the optimal strategy. For a sake of clarity, we show detailed results in table 2 only for three chosen patients. Patient A is in category "very low", patient B in category "low" and patient C in category "very low". We note that for these 3 patients the two cost values are very similar, meaning that we make a good approximation of the value function with our numerical method. We make the same observation on the 47 other patients (data not shown). Also of note is the hierarchy of the cost between the categories of patients. Very low patients have higher optimal costs (between 8.4 and 12) than low (between 3.9 and 9.4) and not too low (between 2.1 and 4.2). This is consistent with the fact that the lower baseline CD4 levels the patient has, the more time will be spent under 500 cells/ µL and the more injections are needed, which both increase the cost of the strategy of injections. We also realized comparisons of several protocols. We simulated five "naive" protocols : P1 with 3-injections cycles, P2 with a first cycle of 3 injections then 2-injections cycles, P3 with 2-injections cycles, P4 with a first cycle of 2 injections then 1-injection cycles and P5 with 1-injection cycles, all protocols with dose 20. Assessing the cost of these protocols is interesting as they imply variable trajectories within the same patient as well as different values for clinical criteria. Moreover they would be clinically feasible and represent a good basis for comparison for our optimal strategy. For every protocol k, we note P +k the space of patients such that cost of optimal strategy is lower than cost of protocol k and n +k its size. We have computed the mean relative positive variation of cost value (MRC), as shown in table 3. We note C opti the mean cost of optimal strategy for patient i and C P ki the mean cost of protocol k for patient i. The MRC allows computing the mean percentage of gain in term of cost function when using the optimal strategy over protocol k:

M RC k = 1 n +k i∈P +k (C P ki -C opti ) C P ki (9)
Results show that mean cost of the optimal strategy is always lower than all other simulated strategies (except one patient for protocol 4, but this is due to numerical approximation, as for this patient W (x 0 )=4.0, C opt = 4.1 and C P4 =4.0). The percentage of cost reduction of the optimal strategy compared to the other protocols in the simulated population of pseudo patients varies from 5.7 to 43%. It confirms that our numerical method allows optimizing the cost function.

In addition to comparing the cost value of all five protocols to the optimal strategy, we have also compared clinical criteria such as the mean number of CD4+ T cells count until horizon, the mean time spent under 500 cells/µL (in days) and the mean number of injections over all simulations. Results of these comparisons are shown in figure 5, where each point corresponds to the value of the criterion for one pseudo-patient, and each color corresponds to the category of the patients ("very low", "low" and "not too low" baseline). We observe that mean cost of the optimal strategy is lower than other simulated protocols and the optimal strategy achieves a good balance between all clinical criteria. Even if CD4+ T cells levels are not as high as for protocols P1, P2 and P3, the optimal strategy allows to spend as much time with levels over 500 cells/µL as these protocols by using less injections. Protocol P5 has the same performance as the optimal strategy for "not too low" patients, as these strategies are very often the same on these patients. The same observation is made on protocol P4 and the "low" patients. Overall, figure 5 shows that the determined strategy allows optimization of the cost function through the chosen criteria (time spent under 500 and number of injections).

More detailed results of comparison of cost function and clinical criteria between optimal strategy and protocols Pk are displayed in table 4 for patients A, B and C. For these three patients, we observe again that mean cost of the optimal strategy is lower than all other simulated strategies. For patient A, the optimal strategy is achieved by two first cycles of two injections then cycles of one injection. For patient B, the optimal strategy consists in a first cycle of 2 injections followed by 1-injection cycles. We can see that the minimum cost is the same for the optimal strategy and protocol P3 (= 5.91): when the patient has a good response to all injections, these strategies are the same. For patient C, the optimal strategy is obtained with 1-injection cycles. Similarly, the minimum cost is the same for the optimal strategy and protocol P5 (= 2.79). For all patients, the optimal strategy is very intuitive: the first complete cycles are needed to raise the number of CD4 over 500 cells/ µL; then, 1injection cycles allow to sustain the levels over 500 cells/ µL. For "not too low" patients, CD4 levels are high enough to use only one injection in the first cycle. This helps reducing the number of injections: in patient A, the optimal Fig. 2 Comparison of cost (3a) and clinical criteria (3b : mean number of CD4+ T cells count until horizon, 3c : mean time spent under 500 cells/µL in days, 3d : mean number of injections over all simulations) between the determined optimal strategy and the five other protocols. Each point corresponds to the value of a pseudo patient, with "very low" patients in yellow, "low" in purple and "not too low" in blue. Mean values within each category are represented by horizontal colored lines. P1: cycles of 3 injections. P2: first cycle of 3 injections then cycles of 2 injections. P3: cycles of 2 injections. P4: first cycle of 2 injections then cycles of 1 injection. P5: cycles of 1 injection. strategy requires 2 to 7 less injections than P1, P2 and P3 but allows to spend as much time over 500 cells/ µL. In patient C the optimal strategy requires one less injection as P3 but allows to spend as much time over 500 cells/ µL. It can be noted that a third injection is never used, even for the first cycles of very low patient. It is due to our choice of cost function: it balances the number of injections and the number of months spent under 500 cells/ µL. The effect of a third injection is usually too low to allow increasing the time spent over 500 cells/ µL by one month and is then not chosen as part of the optimal strategy. These results suggest that our numerical method allows to determine an optimal strategy of injections, and the clinical interpretation of the results are consistent with the mathematical method.

In term of trajectories of the process, figures 3a and 3b show some trajectories obtained with respectively the 2-injection cycles protocol (P3) and the optimal strategy for patient A. We can note that even if CD4 + levels are globally lower in the optimal strategy compared to the two-injections cycles at dose 20 µg/kg, it still allows a maintenance over the threshold of 500 cells/µL by using less injections: indeed, in the best case scenario, the 2-injections cycle strategy implies 5 cycles of 2 injections which is a total of 10 injections, while the optimal strategy induces 2 cycles of 2 injections and 4 single injections, which is a total of 8 injections. The trajectories for patients B and C are provided in appendix D. All together, our results support the interest of determining the optimal strategy based on a criterion combining both the number of injections and the time spent under 500.

Discussion

In this work, we have developed a numerical tool allowing to solve an impulse control problem for a PDMP. The specificity of our work is in the development of a dynamic programming method in the context of a specific biological framework. The objective is to determine the optimal strategy of IL-7 injections for a given HIV-infected patient, in order to maintain CD4 + T lymphocytes levels over the threshold of 500 cells/µL. We first modeled the dynamics of CD4 + T lymphocytes during repeated cycles of IL-7 injections by a PDMP. Then, we solved the impulse control problem by iterating a sequence defined by an integro-differential operator. Theoretical results have shown that this sequence converges to the value function, which allows to determine the optimal action that should be realized at every point of the boundary. We proposed a numerical tool approximating the sequence and the value function on a grid of the state space and applied it to our clinical question. As our method relies on numerical approximation, the obtained optimal strategy could be an approximation of the theoretical one. However the obtained results suggest that we managed to determine optimal strategies for pseudo-patients and that our method allows improving the strategy of injections. Although the horizon of study is only one year, these results are also consistent with a clinical interpretation. The optimal strategy determined for different patients is indeed intuitive: the first cycles aim at increasing the CD4 + T lymphocytes levels and should contain as many injections as possible until the levels are acceptable. Then, the following cycles sustain the CD4 levels over the threshold, and punctual injections are sufficient to reach this objective. The optimal strategy, determined with our method, has a lower cost than other possible clinical strategies. Actually, the obtained optimal strategy depends on the cost previously defined, and we could explore other optimal strategies depending on other cost functions. For example, it could be interesting to use different weights on the time spent under 500 cells /µL and the number of injections (depending on the clinician priorities), or to account for the possible negative side effects due to higher doses (this would need additional data on the question). Finally, the model could be extended by studying the patient until a longer horizon (up to two years). This rises the issue of the increase of computational time (by increasing the size of the grid of the state space) and constitutes a new challenge in itself. In the end, we hope to use this tool in future possible clinical trial investigating the effect of IL-7 injections with patients-specific schedules of injections, personalized and optimized using this method.

applying the results from theorem 1 to determine the optimal cost and an optimal strategy by dynamic programming, we need to determine how to compute numerically the function B to iterate the sequence {Wq} q∈N defined in equation 4. As a reminder, B is defined in [START_REF] Costa | Constrained and unconstrained optimal discounted control of piecewise deterministic markov processes[END_REF] by : K+α)t RV (φ(y, t))dt + e -(K+α)t * (y) TV (φ(y, t * (y)))

BV (y) = [0,t * (y)[ e -(
We will first detail the computation of R then T, and we will finally show how to compute B.

Computation of R

For x = (γ, n, σ, θ, p, r) ∈ X, and function V : X → R, R is defined as :

RV (x) = C g (x) + qV (x) + ηV (x)
with q computing the difference between the states before and after the spontaneous jump.

As Q depends on the action only when the process hits the active boundary,

q(dy|x, d) = η(x)[Q(dy|x) -δx(dy)] = 1 {γ>1} η[δ (1,n,σ,θ,p,r) (dy) -δ (γ,n,σ,θ,p,r) (dy)]
then for every function V , and as K = η :

qV (x) = V (y)q(dy|x) = 1 {γ>1} K[V (1, n, σ, θ, p, r) -V (γ, n, σ, θ, p, r)] Then RV (x) = 1 30 1 {p+r≤500} + qV (x) + KV (x) = 1 30 1 {p+r≤500} + KV (1, n, σ, θ, p, r)1 {γ>1} + KV (x)1 {γ=1} Finally, RV (x) = 1 30 1 {p+r≤500} + KV (1, n, σ, θ, p, r) RV (∆) = KV (∆) (10)

Computation of T

For x ∈ Ξ, and function V : X → R, T is defined as :

TV (x) = inf d∈A(x) C i (x, d) + QV (x, d) = inf d∈A(x) 1 x∈Ξ 1 ∪Ξ 4 + 1 d =0 1 x∈Ξ 3 + V (y) δ (γ(d),1,0,1,Pc,Rc) (dy)1 {x∈Ξ 1 } +δ ∆ (dy)1 {x∈Ξ 2 } + δ (γ(d),n+1,0,θ,p,r) (dy)1 {x∈Ξ 3 } + δ (γ(d),1,0,θ,p,r) (dy)1 {x∈Ξ 4 } +δ (1,n,σ,θ,p,r) (dy)1 {x∈Ξ 5 } Finally, TV (x) = inf d∈A(x) [1 + V (γ(d), 1, 0, 1, Pc, Rc)]1 x∈Ξ 1 +[1 d =0 + V (γ(d), n + 1, 0, θ, p, r)]1 x∈Ξ 3 +[1 + V (γ(d), 1, 0, θ, p, r)]1 x∈Ξ 4 + V (∆)1 x∈Ξ 2 + V (1, n, σ, θ, p, r)1 x∈Ξ 5 TV (∆) = V (∆) (11)

Computation of B

Now, for Y ∈ X, and function V : X → R, we need to compute : K+α)t RV (φ(y, t))dt + e -(K+α)t * (y) TV (φ(y, t * (y)))

BV (y) = [0,t * (y)[ e -(
As we cannot make an exact computation of BV on X, we need to approximate this computation on a grid of the state space. In order to detail the approximation of the computation, we define K+α)t RV (φ(y, t))dt and H(V, y) = e -(K+α)t * (y) TV (φ(y, t * (y))

G(V, y) = [0,t * (y)[ e -(
as in 7 and 8. We define a time interval ∆t (in practice equal to one day) and for every y = (γ, n, σ, θ, p, r) ∈ X, we note n * (y) = t * (y) ∆t

For every j ∈ {0..n * (y)-1}, we note φ j (y, t) = φ(y, j∆t) and φ(y, t * (y)) = (γ, n, σ+t * (y), θ+ t * (y), p * (y), r * (y)). The integral defined in equation 7 is computed by approximation using the classic trapezoidal rule using the j∆t nodes :

G(V, y) ∆t 2 RV (y) + ∆t 2 e -(K+α)t * (y) RV (φ(y, t * (y))) + n * (y)-2 j=1 ∆te -(K+α)j∆t RV (φ j (y, t))
with RV (x) = 1 30 1 {p+r≤500} + KV (1, n, σ, θ, p, r), as computed in 10. Then we obtain the following for every y = (γ, n, σ, θ, p, r) ∈ X : e -(K+α)j∆t ( 1 30 1 {p j +r j <500} + KV (1, n, σ + j∆t, θ + j∆t, p j , r j ))

G(V, y) = ∆t 2 ( 1 30
1 {p+r<500} + KV (1, n, σ, θ, p, r)) + ∆t 2 e -(
Now, we need to compute H as defined in 8 : it depends on TV (φ(y, t * (y))), which takes different values according to the boundary reached in that point, as written in equation 11.

Moreover, as we know the flow, we can give conditions on y = (γ, n, σ, θ, p, r) to reach a given boundary in φ(y, t * (y)). Then :

-if φ(y, t * (y)) ∈ Ξ 1 (θ ≤ 1) then H(V, y) = inf d∈[d 1 ,..dm d ] e -(K+α)t * (y) 1 + V (γ(d), 1, 0, 1, Pc, Rc) -if φ(y, t * (y)) ∈ Ξ 2 (θ + t * (y) ≥ T h ) then H(V, y) = e -(K+α)t * (y) V (∆) -if φ(y, t * (y)) ∈ Ξ 3 (n < n inj , θ + t * (y) < T h ) then H(V, y) = inf d∈[0,d 1 ,..dm d ] e -(K+α)t * (y) 1 {d =0} +V (γ(d), n + 1, 0, θ + t * (y), p * (y), r * (y)) -if φ(y, t * (y)) ∈ Ξ 4 (n = n inj , γ = 1, θ + t * (y) < T h ) then H(V, y) = inf d∈[d 1 ,..dm d ]
e -(K+α)t * (y) 1 + V (γ(d), 1, 0, θ + t * (y), p * (y), r * (y))

-if φ(y, t * (y)) ∈ Ξ 5 (n = n inj , γ > 1, θ + t * (y) < T h ) then H(V, y) = e -(K+α)t * (y) V (1, n, σ + t * (y), θ + t * (y), p * (y), r * (y))
Finally, for every y = (γ, n, σ, θ, p, r) ∈ X : The grid of the state space created in Matlab contains 67614 lines and 7755 columns. For a given patient, the computation of 40 iterations of the sequence (convergence is reached between 35 and 45 iterations) requires between 5 and 6 days.

BV (y) = ∆t 2 ( 1 30
1 {p+r<500} + KV (1, n, σ, θ, p, r)) + ∆t 2 e -(K+α)t * (y) 1 30 1 {p * +r * <500} + KV (1, n, σ + t * , θ + t * , p * (y), r * (y)) +∆t n * (y)-2 j=1 e -(K+α)j∆t 1 30 1 {p j +r j <500} + KV (1, n, σ + j∆t, θ + j∆t, p j , r j ) + inf d∈[d 1 ,..dm d ] e -(K+α)t * (y) 1 + V (γ(d), 1, 0, 1, Pc, Rc) 1 {θ≤1} +e -(K+α)t * (y) V (∆)1 {θ+t * (y)≥T h } + inf d∈[0,d 1 ,..dm d ] e -(K+α)t * (y) 1 {d =0} +V (γ(d), n + 1, 0, θ + t * (y), p * (y), r * (y)) 1 {n<n inj ,θ+t * (y)<T h } + inf d∈[d 1 ,..dm d ] e -(K+α)t * (y) 1 +V (γ(d), 1, 0, θ + t * (y), p * (y), r * (y)) 1 {n=n inj ,γ=1,θ+t * (y)<T h } +e -(K+α)t * (y) V (1, n, σ + t * (y), θ + t * (y), p * (y), r * (y))1 {n=n inj ,γ>1,θ+t * (y)<T h } (12) and BV (∆) = [0,+∞) e -(K+α)t KV (∆)dt = K K + α V (

C Sensitivity analysis of the method

To evaluate how the uncertainty on individual parameters estimation could impact the determination of the optimal strategy, we have realized a sensitivity analysis. For a given patient, we suppose a normal distribution of parameters λ and ρ. We generate L = 500 pairs of parameters (λ, ρ) from this joint distribution. Each pair corresponds to an initial

}

Step 1: Determine x vs = (γ, n, σ, θ, p, r) = χ -1 ((v, s)) by identifying block i (γ,n), line j (σ,θ) using S1, column s (p,r).

Step 2: Get t * (x vs ) with S3 and the flow from x vs with S2 =⇒ define the elements of the grid x m necessary for the computation of GW q (x vs ).

Step 3: Determine the positions χ(x m ) of x m in the grid using S1.

Step 4: From previous iteration, get W q (x m ) = M q (χ(x m )) and use linear combination determined in Appendix A to compute GW q (x vs )

Step 5: Determine boundary reached in t * (x vs ) and compute HW q (x vs ) as in value of lymphocytes T CD4 0 . We determine the empirical quartiles of the distribution of the CD4 0 and focus on the pairs inducing values close the first and the third quartiles.

Then, for each pair we simulate the five possible protocols P1 to P5 and compare them to the optimal strategy determined on the mean value of (λ,ρ). In practice, values of pairs and associated values of CD4 are displayed in table 5. For the mean value of (λ,ρ), we determined the optimal strategy to be a first cycle of 2 injections and then cycles of 1 injection, which corresponds to protocol P4. We show in table 5 the cost of each protocol for each pair of (λ,ρ), and we put in bold the minimum cost over the five protocols. We can see that protocol P4 achieves the minimum cost for all pairs inducing CD4 values at the first quartile. For pairs inducing CD4 values at the third quartile, the protocol achieving the minimum cost is P5. However, the difference of cost is not huge and P4 actually induces more time spent over the 500 threshold and less than one more injection than P5 on average, which is still acceptable. Overall, this shows that even with some error on the estimation on λ, ρ we would be able to determine a strategy achieving a good balance between clinical criteria.

D Trajectories of patients B and C

28 Chloé Pasin et al.

Fig. 1

 1 Fig. 1 Mechanistic model for the dynamics of CD4 + T lymphocytes.

  Dynamics of CD4 + T lymphocytes in patient A under P3, a 2-injections cycles protocol (dose 20). Dynamics of CD4 + T lymphocytes in patient A under the determined optimal strategy.

Fig. 3

 3 Fig.3Dynamics of CD4 + T lymphocytes in patient A. Straight line corresponds to the "best" outcome, i.e., when the effect of all injections lasts seven days. Dashed line corresponds to other possible trajectories, when this effect can last less than seven days.

  ∆) B Structure of the code Structure of the code and its subroutines are shown in figure 4. Application in the results section requires the following grid : γ ∈ {1..3} n ∈ {1..3} σ ∈ {0..351} θ ∈ {1..365} p ∈ {2..110} depending on the patient r ∈ {100..1500} depending on the patient

  Dynamics of CD4 + T lymphocytes in patient B under P3, a 2-injections cycles protocol (dose 20). Dynamics of CD4 + T lymphocytes in patient B under the determined optimal strategy.

Fig. 5

 5 Fig.5Dynamics of CD4 + T lymphocytes in patient B. Straight line corresponds to the "best" outcome, i.e., when the effect of all injections lasts seven days. Dashed line corresponds to other possible trajectories, when this effect can last less than seven days.

Fig. 6

 6 Fig.6Dynamics of CD4 + T lymphocytes in patient C. Straight line corresponds to the "best" outcome, i.e., when the effect of all injections lasts seven days. Dashed line corresponds to other possible trajectories, when this effect can last less than seven days.

Table 1

 1 Characteristics of the pseudo-patients population.

	Parameter	Mean (std)
	λ (cells/day) 2.24 (0.39)
	ρ (/day)	1.96 (0.84)
	π 0 (/day)	0.0461 (0.0035)
	µ R (/day)	0.0503 (0.0033)
	µ P (/day)	0.0717 (0.014)
	βπ 1	0.958
	βπ 2	0.752
	βπ 3	0.143
	Category	Number of patients (%)
	Very low	4 (8%)
	Low	24 (48%)
	Not too low	22 (44%)

Table 2

 2 Comparison of cost values from value function and Monte Carlo simulation. Std = standard deviation.

		Patient A	Patient B	Patient C
		very low	low	not too low
	Optimal cost W (x 0 )	9.47	6.11	2.87
	Cost of optimal strategy : mean (std) 9.53 (0.85) 6.20 (0.56) 2.90 (0.36)
	(obtained by Monte Carlo)

Table 3

 3 Computation of the mean relative variation of cost value (MRC) for every protocol allows determining the mean percentage of gain in term of cost function when using the optimal strategy over protocol k. P1: cycles of 3 injections. P2: first cycle of 3 injections then cycles of 2 injections. P3: cycles of 2 injections. P4: first cycle of 2 injections then cycles of 1 injection. P5: cycles of 1 injection.

	Protocol	P1	P2	P3	P4	P5
	n +	50	50	50	49	50
	MRC	43%	31% 20% 5.7% 8.8%

Table 4

 4 Comparison of protocols of injections for patient A, B and C. P1: cycles of 3 injections. P2: first cycle of 3 injections then cycles of 2 injections. P3: cycles of 2 injections. P4: first cycle of 2 injections then cycles of 1 injection. P5: cycles of 1 injection.

	Protocol	P1	P2	P3	P4	P5	Optimal
	Patient A
	Mean cost	14.8 11.6 10.7 10.4 11.1	9.53
	Std	0.52 0.62 0.68 1.21 1.23	0.85
	Min cost	14.6 11.4 10.4 9.56 10.1	8.82
	CD4 mean	671	662	659	552	506	578
	Days under 500	54.5 55.4 56.5 102	150	58.0
	Number of injections 15.1 11.1 10.2 8.25 7.48	8.66
	Patient B
	Mean cost	9.35 7.62 6.91 6.26 6.63	6.20
	Std	1.14 1.00 0.98 0.70 0.87	0.56
	Min cost	8.62 6.91 6.14 5.91 6.12	5.91
	CD4 mean	762	742	736	622	598	625
	Days under 500	18.4 19.0 24.2 28.1 61.9	24.4
	Number of injections 9.83 7.73 6.84 6.05 5.27	6.10
	Patient C
	Mean cost	5.54 4.75 3.80 3.08 2.93	2.90
	Std	0.29 0.24 0.31 0.37 0.49	0.36
	Min cost	5.46 4.70 3.73 2.95 2.79	2.79
	CD4 mean	774	762	758	666	631	631
	Days under 500	5.56 5.53 5.65 5.90 8.43	5.89
	Number of injections 6.02 5.02 4.04 3.13 3.03	3.07

Table 5

 5 Pairs of (λ,ρ), associated CD4 values and mean cost for protocols P1 to P5. P1: cycles of 3 injections. P2: first cycle of 3 injections then cycles of 2 injections. P3: cycles of 2 injections. P4: first cycle of 2 injections then cycles of 1 injection. P5: cycles of 1 injection. P4 is the optimal protocol for mean value of (λ,ρ).

	Category λ	ρ	CD4 0	P1	P2	P3	P4	P5
	Mean	2.065 2.022 289	6.55 5.54 5.54 4.16	4.95
	Q1	1.506 2.305 224	8.59 6.94 6.03 5.20	5.67
	Q1	2.062 1.180 223	11.6 9.19 8.52 7.28	8.02
	Q1	1.701 1.747 224	8.94 7.31 6.50 5.60	6.40
	Q1	2.163 1.078 222	11.9 9.50 8.82 8.04	8.39
	Q1	1.737 1.371 224	9.06 7.39 6.64 5.71	6.48
	Q1	1.689 1.758 223	8.92 7.31 6.52 5.62	6.40
	Q1	1.728 1.689 224	9.00 7.35 6.63 5.71	6.45
	Q1	1.426 2.599 222	8.54 6.90 5.98 5.14	5.62
	Q1	2.493 0.838 222	13.3 11.0 10.3 8.85	9.76
	Q1	1.805 1.542 223	9.23 7.59 6.94 6.36	6.66
	Q3	2.160 2.594 337	5.62 4.81 3.88 3.15	3.04
	Q3	2.424 1.956 336	5.83 4.99 4.10 3.81	3.21
	Q3	2.638 1.625 335	6.10 5.22 4.30 3.93	3.31
	Q3	2.219 2.429 337	5.68 4.85 3.93 3.21	3.11
	Q3	2.477 1.879 337	5.87 5.04 4.12 3.83	3.21
	Q3	2.466 1.896 337	5.86 5.01 4.09 3.82	3.22

Acknowledgements We would like to thank the main investigators and supervisors of INSPIRE 2 and 3 studies: Jean-Pierre Routy, Irini Sereti, Margaret Fischl, Prudence Ive, Roberto F. Speck, Gianpiero D'Ozi, Salvatore Casari, Sharne Foulkes, Ven Natarajan, Guiseppe Tambussi, Michael M. Lederman, Therese Croughs and Jean-François Delfraissy. This work was supported by the Investissements d'Avenir program managed by the ANR under reference ANR-10-LABX-77.

Appendix A.

Step 6:

Compute W q+1 (x vs ) = GW q (x vs ) + HW q (x vs )

Outputs: at each iteation q, a matrix M q of size N sum × N pr with value of W q (x vs ) in position (v, s).

Subroutine 1 (S1)

Creation of combinations (σ,θ) Outputs: for every combina-

Creation of t * (x) and W 0 (x) for every x of the grid, identified using S1.

Outputs: a matrix of size Fig. 4 Structure of the code and its subroutines.


Nermine Chaari

David Rodney

Emmanuel Clouet
email: emmanuel.clouet@cea.fr

Oxygen -dislocation interaction in titanium from first principles

Keywords: Dislocation, Plasticity, Titanium, Oxygen, Hardening

come

Oxygen-dislocation interaction in titanium from first principles

In hexagonal close-packed (hcp) titanium alloys (α-Ti alloys), plasticity is strongly influenced by the interactions between 1/3 1210 dislocations and interstitial solute atoms, in particular oxygen [1][2][3][4][START_REF] Caillard | Thermally activated mechanisms in crystal plasticity[END_REF]. A small oxygen content leads to an important hardening, with the plastic flow controlled by thermal activation instead of the almost athermal plasticity observed in pure titanium. Such hardening cannot be explained by a simple elastic interaction between dislocations and impurities [1,4,[START_REF] Yu | [END_REF]. The most likely reason is a change of dislocation core properties because microscopic observations and in situ experiments show long rectilinear screw dislocations in titanium with oxygen [4,[START_REF] Yu | [END_REF][7][8], an evidence that the lattice friction against the glide of screw dislocations increases in presence of oxygen atoms and may control the plasticity.

Transmission electron microscopy (TEM) also shows that dislocation motion is intermittent at low temperature, with long time periods where the lines are immobile, followed by short events where the dislocations glide in prismatic planes over a long distance [7,9,10]. This motion becomes more jerky when the oxygen content increases, with a decrease of the glide distance in prismatic planes [8].

Dislocations also cross-slip at higher rates, mostly from prismatic to first-order pyramidal planes, but also to basal planes for well-oriented samples [START_REF] Churchman | Proc. R. Soc. Lond. A[END_REF]. Jogs and super-jogs appear on screw dislocations because of multiple cross-slip, finally leading to the creation of dislocation dipoles and debris [START_REF] Yu | [END_REF][7][8][START_REF] Naka | [END_REF][13][14].

High-resolution electron microscopy gives some information on the dislocation core structure in titanium alloys [START_REF] Yu | [END_REF]15,16]. In particular, a narrowing of the screw core has been evidenced in presence of oxygen [START_REF] Yu | [END_REF]. Ab initio calculations based on the Density Functional Theory have also proved a suitable approach [17]. It was shown that the ground state of the screw dislocation in pure titanium is dissociated in a first-order pyramidal plane and needs to transform into a metastable configuration dissociated in a prismatic plane to glide easily [10], leading to a locking-unlocking mechanism, in agreement with the intermittent glide observed in TEM in situ straining experiments [7,9,10]. Moreover, previous calculations [START_REF] Yu | [END_REF] have evidenced a strong short-range repulsion between the prismatic core of the screw dislocation and an oxygen atom. Here, we extend this work to consider all possible configurations of the screw dislocation core and in particular the pyramidal ground state, in order to determine how oxygen atoms affect transitions between core configurations and what are the potential impacts on the plasticity of Ti-O alloys. Ab initio calculations are performed with the Pwscf code [18], in the generalized gradient ap-proximation using the PBE functional [19]. Ultrasoft pseudopotentials of Vanderbilt type are used with 3s and 3p semi-core electrons included in valence states for titanium. Electronic wavefunctions are expanded in plane waves with a cutoff energy of 28 Ry. A regular grid of 2 × 1 × 7 k-points is used for integration in reciprocal space and the electronic density of states is broadened with the Methfessel-Paxton function using a spreading of 22 mRy. Atoms are relaxed until all atomic force components are below 10 meV/ Å.

To preserve periodic boundary conditions, we introduce a dislocation dipole in the simulation cell with a quadrupolar geometry [17,20]. The periodicity vectors, before introduction of the dipole, are a x = 6a[1010], a y = 8c[0001], and a z = 3 b = a[1210], corresponding to 576 lattice sites. This setup correctly predicts that the dislocation ground state is dissociated in a first-order pyramidal plane [10] (Fig. 1a). The prismatic core (Fig. 1f) has a positive relative energy ∆E = 5.96 meV/ Å. We also consider another metastable configuration (Fig. 1j) [10], which is spread in both prismatic and pyramidal planes and has a higher relative energy, ∆E = 11.0 meV/ Å. For the oxygen atom, our calculations predict that the most stable insertion site is octahedral (O), with excess energies for the other metastable positions in agreement with previous ab initio studies [21] (see Supplementary material).

To study the interaction between screw dislocations and oxygen atoms, we start from a dislocation relaxed in pure titanium in one of the three reference configurations. We then an oxygen atom in the vicinity of the dislocation core. Since the simulation cell is three Burgers vectors long in the z direction, we effectively insert an oxygen atom every three sites along the dislocation line. The interaction energy is defined as

E int = E dislo-O -E O -E dislo + E bulk , (1)
where E dislo-O , E O , E dislo and E bulk are the energies of the same supercell with respectively both the dislocation dipole and the oxygen atom, with only the oxygen atom, only the dislocation dipole, and no defect nor solute. Attraction corresponds to a negative interaction energy. The dislocation reference energy E dislo is that of the initial dislocation configuration before insertion of the oxygen atom. As only one dislocation is decorated with an oxygen atom, part of this interaction energy arises from a variation of the elastic interaction between the two dislocations when the interaction with the oxygen atom is strong enough to induce a motion of the dislocation. Relaxed configurations are then analyzed with the help of differential displacement maps [22] and the screw component of the Nye tensor [23]. The results of our calculations are summarized in Fig. 1, where each line corresponds to one of the three initial core configurations and the columns show the relaxed cores, without oxygen atoms in the first column and with oxygen atoms inserted in different positions in the other columns.

We first consider the case when an oxygen atom is inserted in the dislocation core in an octahedral position of the perfect crystal that is destroyed by the stacking fault ribbon when the dislocation dissociates [START_REF] Yu | [END_REF][24][25][26]. These sites are shown as diamonds in Fig. 1 in the pyramidal and prismatic planes. These positions are strongly repulsive for all three initial core configurations, which all become unstable and relax to new dislocation cores. The ground state, initially dissociated in a pyramidal plane (Fig. 1a), relaxes to a compact core (Fig. 1b,c). The same compact core is obtained with the metastable prismatic configuration (Fig. 1f) when the oxygen atom is inserted at the dislocation center (Fig. 1h). For an insertion site away from the dislocation center, the prismatic configuration is destabilized and falls back into the pyramidal ground state (Fig. 1g). The same destabilization is obtained with the mixed pyramidalprismatic metastable core (Fig. 1j), which also reverts back to the pyramidal core (Fig. 1k). In that case, the interaction energy is negative (E int = -80 meV) but this is mainly due to the energy variation associated with the change of core configuration (3b ∆E = 97 meV).

The compact core identified in Fig. 1b, c and h is of high energy and can therefore be produced only under an applied stress. A narrowing of the pyramidal and prismatic spreading induced by oxygen atoms however agrees with the high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) observations of Yu et al. [START_REF] Yu | [END_REF], who reported that the screw dislocation core is smaller in a titanium alloy containing 0.3 than 0.1 wt. % oxygen, with oxygen atoms occupying octahedral sites in the core vicinity. Similar highenergy compact cores are also induced by substitutional solutes, in particular indium [27].

The second case considered here is when the oxygen atom is introduced in insertion sites created by the stacking fault ribbon when the dislocation dissociates. Interstitial sites O P (square symbols in Fig. 1) can be found in both the prismatic and pyramidal faults. Our calculations show that these insertion sites are stable and do not induce any core reconfiguration (Fig. 1d,i,l). The O P sites nevertheless lead to a repulsive interaction, with interaction energies in the range 100 -140 meV. These interaction energies are higher than reported by Yu et al. (50 meV) [START_REF] Yu | [END_REF] but the difference may be because these authors used a cell with a smaller height h = b, thus corresponding to the dislocation interaction with a full row of oxygen atoms.

The third case is when the stacking fault ribbon lies in the pyramidal plane and the oxygen atom is inserted in a new O π octahedral site created by the pyramidal fault (triangle symbols in Fig. 1). These interstitial sites correspond to octahedral sites in the twined hcp crystal [26], since the pyramidal fault is a two-layer twin in the pyramidal stacking [28,29]. Both the pyramidal and mixed pyramidal-prismatic cores relax towards the same core shown in Figs. 1e andm, which is close to the mixed pyramidal-prismatic core in pure titanium with the oxygen atom at the center of the core. Compared to the pyramidal core, the interaction is again repulsive (E int = 87 meV).

From the above calculations, we conclude that the lowest-energy configuration in Ti-O alloys remains the pyramidal core without oxygen atoms. All other configurations, although metastable, have a higher energy. Therefore, we can anticipate that in Ti-O alloys, oxygen atoms do not segregate on screw dislocations. To glide, the pyramidal core needs to transform into a prismatic core. When the latter comes across an oxygen atom in an octahedral site, i.e. a site that will be destroyed by the fault, the first site encountered by the dislocation leads to the pyramidal core in Fig. 1g. The compact core in Fig. 1h is less likely because of its much higher energy and because it corresponds to an octahedral site further inside the dislocation core. The decorated prismatic core in Fig. 1i is also unlikely, not only because of its high energy, but also because it requires migration of the oxygen atom to a newly created octahedral site while the dislocation is in an excited state and glides in the prismatic plane.

At low temperatures, when oxygen long-range diffusion is not active, the most likely scenario is therefore that, when the glissile prismatic core encounters an oxygen atom, it locally reverts back to a pyramidal core to avoid the oxygen atom. This scenario leads to dislocation cross-slip and is thus specific to the screw orientation, in agreement with experiments [4,[START_REF] Yu | [END_REF][7][8], which report that screw dis-locations are more affected by the presence of oxygen atoms than edge dislocations. This scenario is also consistent with TEM observations [8], which show that an oxygen addition destabilizes prismatic slip, leading to shorter glide distances in prismatic planes between pinning points. As the gliding dislocation locally cross-slips to bypass oxygen atoms, jogs are created [START_REF] Yu | [END_REF]26] in agreement with the postmortem TEM observations of Williams et al. [30] showing that in the dislocation microstructure that develops in presence of oxygen atoms, most dislocations are jogged screw dislocations. These jogs are responsible for an increased lattice friction acting against dislocation glide. They also lead in αtitanium alloys to creep rates controlled by the motion of jogged screw dislocations [31,32]. When the temperature becomes high enough to allow for oxygen diffusion, this bypass mechanism will be in competition with oxygen migrating out of the repulsive region, thus allowing the dislocation to remain in its glide plane.

A reconstruction of the dislocation core into a pyramidal core induced by the oxygen atoms also allows to understand why cross-slip in pyramidal planes becomes more active with an increase of the oxygen content [START_REF] Yu | [END_REF]8]. Indeed, once stabilized, the pyramidal core can glide in its pyramidal habit plane by nucleating kink pairs [10].

Comparing the present results with Ref. [26] on the interaction of screw dislocations with oxygen atoms in zirconium, we see that the same repulsive interaction resulting from the destruction of oxygen insertion sites by the stacking fault ribbon and leading to a change of dislocation core configuration is obtained in zirconium and titanium when the oxygen atom lies in its natural octahedral insertion sites. Therefore, in both metals, the screw dislocations gliding in prismatic planes undergo a stronger lattice friction because of the creation of jogs allowing to bypass the solute atoms. Moreover, cross-slip in pyramidal planes is enhanced in both metals by oxygen atoms. However, while in titanium, cross-slip in pyramidal planes is already active at low temperatures even for a very low oxygen content [10,33], in zirconium, plasticity is confined to prismatic planes at low temperatures and crossslip appears only at higher temperatures thanks to thermal activation or interaction with oxygen atoms [26]. The origin of the lattice friction in the principal slip system also differs between both metals: in zirconium, this friction is purely extrinsic and only arises from the interaction with impurities, in particular oxygen, whereas an intrinsic friction still remains in pure titanium at low temperature [10,33,34].

Figure 1 :

 1 Figure 1: Core structure of a dislocation initially dissociated (a-e) in the pyramidal plane, (f-i) in the prismatic plane, and (j-m) in both the pyramidal and prismatic planes, in presence of an oxygen impurity. The initial configurations of the dislocation, before introduction of the oxygen atom, are shown in (a, f, j) with the various insertion sites for the oxygen atom shown in green. Relaxed structures are shown in (b, c, g, h, k) for an oxygen atom lying in an octahedral site destroyed by the stacking fault ribbon of the dislocation, in (d, i, l) for an oxygen atom in an octahedral site created by the prismatic component of fault, and in (e, m) by its pyramidal component. In these projections perpendicular to the dislocation line, titanium atoms are sketched by symbols with a color depending on their (1210) plane in the original perfect crystal and with different symbols used for atoms, depending on their neighborhood in the dislocated crystal. The arrows between two atomic columns are proportional to the differential displacement created by the dislocation in the [1210] direction. Displacement smaller than 0.1 b are not shown. The contour map shows the dislocation density according to the Nye tensor.

Acknowledgements -This work was performed using HPC resources from GENCI-CINES and -TGCC (Grants 2017-096847). The authors also acknowledge PRACE for access to the Curie resources based in France at TGCC (project PlasTitZir). DR acknowledges support from LABEX iMUST (ANR-10-LABX-0064) of Université de Lyon (program "Investissements d'Avenir", ANR-11-IDEX-0007)


Armand Le Gouguec
email: armand.le-gouguec@renault.com

Andras Kemeny
email: andras.kemeny@renault.com

Alain Berthoz
email: alain.berthoz@college-de-france.fr

Frédéric Merienne
email: frederic.merienne@ensam.eu

Artificial Potential Field Simulation Framework for Semi-Autonomous Car Conception

Keywords: Artificial potential field, Autonomous car, Simulation, Human-machine interface

Artificial potential field is investigated to provide a high level of synergy between driver and semiautonomous vehicle. This article presents a framework developed to test the performances of this approach. Stand-alone performances of this system is tested for a lane keeping and cruise control application. Performances are promising and future development is discussed.

Introduction

Autonomous cars are one of the next major innovations in the field of car manufacturing. Such technology represents a real revolution, especially regarding safety, labor organization and town planning [Eug13]. However, driving is a difficult task to automatize due to its highly dynamic nature and the variety of situation in which it can be exercised. Therefore, most of the time, autonomous driving is approached in a progressive way. However, autonomous vehicles initially designed to operate under specific condition can rarely be adapted to operate in less specific situations. Also, autonomous vehicles only handling some situations require an autonomous-to-manual transition so a human driver can take over controls when leaving the autonomous system comfort zone. And partially autonomous vehicles able to execute only some specific part of the driving task need to be assisted by a human operator at all time. In both cases, a good cooperation between driver and autonomous car is required. Whether it is to share vehicle control or to switch between autonomous and manual driving, the driver needs to understand the intentions of the autonomous system quickly and correctly. In pursuit of a solid base for polyvalent autonomous vehicle development, as well as for efficient humancomputer interaction, this work explores the use of artificial potential field method applied to autonomous driving technology by proposing an architecture for autonomous vehicle low-level control.

Theoretical background

Artificial Potential Field (APF) is a method used to make robots navigate in an environment containing obstacles [Kha86,Oku86]. Originally, this method was developed to control robotic arms in real time by means of force fields. Those fields were designed to attract the robot into the desired place and to keep it away from any obstacles at the same time as seen Fig. 1. APF method differs from classic motion planning technique due to its reactive nature contrasting with traditionally more proactive methods. For this reason, APF method requires no preliminary computation to run. Therefore, environment modification or dynamic obstacles do not have significant impact on performances which gives the APF method an edge over classic motion planning for highly dynamic environments. Classic APF guidance is based on two types of fields. Attractive field attached to the goal and repulsive field attached to obstacles. Once generated, those fields are summed up together in order to generate the global guidance field used for the robot navigation. The field generated by this operation is then used to apply guidance forces to the robot using gradient descent. The force obtained by this process can be applied as acceleration setpoint to control the robot.

Even if the APF approach is better suited for holonomic robots, it can be adapted to autonomous driving problems. Semsar-Kazerooni applied some of the APF method to create an automated cruise control enhance with platooning capability [SK16]. In this system, longitudinal distance between vehicles was maintained by a specific artificial potential field.

APF technology was also used to produce racing artificial intelligence. By mixing multi agent system with artificial potential field method, Uusitalo managed to build an efficient autonomous racing car system even if purely reactive [Uus11].

Wolf implement an autonomous system able to drive on the highway [Wol08]. This also proves APF method can provide strategic decision-making up to a certain limit.

Research question

The reason behind focusing on APF guidance come from the similarities between APF and human cognition models. This has been theorized by Gibson [Gib38] as shown by Fig. 2 and confirmed by Fajen [Faj03a,Faj03b]. This parallel between cognition and guidance fields makes APF a good candidate to realize a semi-autonomous vehicle easy to empathize with.

In order to experiment with this hypothesis, a test environment is required. Our question of research is related to the creation of a simulation framework able to handle complex dynamic model and various field representation of the different obstacles. This framework also needs to run real time to enable human-in-the loop testing.

Our proposed architecture integrates the APF method in a way shown in Fig. 3. Important entities are extracted from the environment and translated in terms of potential field. Those fields are then merged together to obtain the final guidance field. On top of the APF generation is built a control loop designed to interact with the car and the driver as shown Fig. 4. This system is responsible for converting forces into commands for the vehicle as well as interacting with the driver.

Development

In order to illustrate the internal functioning of the framework, a simple implementation example is shown in this section. The situation considered here implies a sole autonomous vehicle which aims to maintain speed while staying in the middle of the lane. Implementation is made using the modular driving simulation software SCANeR studio (version 1.6).

In this context, the only obstacle that needs to be considered is the edge of the current lane. Because this is a simulation, all the information about the environment is directly accessible, thus reducing the perception task to extracting the road geometry through the SCANeR API. Fig. 5 show what the information extracted during the perception process can look like. Classic obstacle fields are attached to the edges of lanes in order to push the vehicle back toward the center of the road. Those fields are defined by Eq. 1.

𝑈 𝑜 (𝑋 ⃗ ) = { 1 2 𝜂 ( 1 𝜌(𝑋 ⃗⃗) - 1 𝜌 0 ) 2 𝑖𝑓 𝜌(𝑋 ⃗ ) ≤ 𝜌 0 0 𝑖𝑓 𝜌(𝑋 ⃗ ) > 𝜌 0 (1)
η is a constant which is used to modulate the field intensity. With 𝑋 ⃗ being the robot position in space, ρ(𝑋 ⃗ ) represents the distance between the robot and the obstacle. ρ 0 is the limit radius of the repulsive field area of effect. Beyond this distance, the obstacle doesn't interact with the robot anymore. Details about those fields can be found in [Kha86]. Fusion of the different fields is done by summing up all fields obtained in the previous step.

In order to calculate the force required, the guiding field is then sampled at multiple locations under the vehicles and inclination is determined through planar least square regression. An additional artificial friction force is also introduced to limit oscillations. This force is also responsible for cruise control. Artificial friction is defined by Eq. 2.

𝐹 ⃗ 𝑓 (𝑉 ⃗⃗ ) = 𝜇(𝑉 ⃗⃗ 0 -𝑉 ⃗⃗ )

(2)

V ⃗ ⃗⃗ is the current velocity of the vehicle. V ⃗ ⃗⃗ 0 is the nominal velocity of the road. This nominal velocity is oriented in the same direction as traffic and its value is the speed required by the driver. µ is a constant used to modulate the intensity of the force. The sum of the artificial friction force with the previously computed lane keeping force gives the global force to apply to the vehicle. Control is achieved through a reverse control maps. Those maps give controls to apply for a given system status and acceleration setpoint. Reverse control maps are generated by benchmarking the vehicle performances prior to the experiment. As an example, throttle reverse map is shown Fig. 6.

Results and Discussion

Tests were conducted on a terrain used to validate some autonomous vehicle capabilities. It consists of a four-lane highway displaying a succession of curves with a radius of curvature of 400m. The road layout is shown Fig. 7. The vehicle behaves as intended for a wide range of parameters. However, extreme setups can lead to some erratic behavior. Insufficient artificial friction leads to system instability resulting in oscillations and eventually lane departures. A small field area of effect leads to violent force variations and unpredictable behavior. The same phenomenon is observed if the repulsion gain is too low.

As an example of "normal" behavior, following results are obtained for a required speed of 100 km/h and with gain parameters η and µ both set to 10. The obstacle effect radius ρ 0 is set to 3 meters and forces are computed by regression over 4 points placed in the corners of a 30 cm square around the vehicle center of gravity.

In this setup, vehicle speed is maintained within 1% of the requested speed and acceleration never rises above 2 m/s². At the corners, the vehicle shifts toward the outside as if centrifugal force was pushing it as shown Fig. 7. This is because the chosen obstacle field profile doesn't take speed into account.

This issue might be reduced or negate by the use of dynamic potential fields adapting shapes according to the relative velocity between the vehicle and the obstacle. An example of a dynamic potential field can be found in [Par08]. A maximum of 40 cm offset was recorded as shown Fig. 8. Considering that lanes are 3.7 m wide, this lateral shifting is somewhat acceptable even if not desired. The framework and its implementation don't create any noticeable perturbation in regard to the APF method. Reverse control maps performed extremely well by keeping acceleration within 0.1% of the setpoint despite being an open loop controller.

However, the dynamic model used was fairly simple and more complex ones might impact control performances. The whole system was running in real time at 100 Hz on a laptop without any performances issue. This frequency was sufficient for the tested situation but more complex environment and haptic feedback may require higher execution speed.

Conclusion and Perspectives

This work lays the foundations of an architecture which aims to provide a new approach for autonomous vehicles conception based on human machine cooperation and artificial potential fields. Stand-alone performances in simple situations are promising and more exhaustive testing are on the way including in depth parametric exploration and testing of other artificial potential field profiles.

Corresponding development of the human-machine interface is also in progress starting with real time graphical representation of the global guiding field.

In addition to being an efficient tool for designing potential fields, it might be an intuitive way to convey the intentions of the autonomous system to the driver in addition to force feedback. If future human experimentations confirm this hypothesis, the architecture described here will take on its full significance as an efficient and user-friendly paradigm.

Figure 1 .

 1 Figure 1. Contour plot of an artificial potential field and trajectory of a robot guided by it. Black polygons represent obstacles. The point named S is the starting position and G is the goal. The trajectory is obtained by a gradient descent procedure along the field. Figure from Okutomi's work [Oku86]

Figure 2 .Figure 3 .

 23 Figure 2. A field of safe travel at a traffic intersection. An attempt to represent how proximity to other vehicles trajectories affect our own driving behavior. Illustration from Gibson's work [Gib38].

Figure 4 .

 4 Figure 4. Illustration of the main interactions between the different parts of the framework.

Figure 5 .

 5 Figure 5. SCANeR graphic layer (a) and road logic extracted (b) side by side.

Figure 6 .

 6 Figure 6. Throttle reverse control map.

Figure 7 .

 7 Figure 7. Road layout used for the test and vehicle trajectory obtained. Lateral shifting is magnified (150x) to show the effects of road curvature on the trajectory.

Figure 8 .

 8 Figure 8. Lateral shifting along the track.