Camille Pouchol
email: pouchol@ljll.math.upmc.fr

Jean Clairambault
email: jean.clairambault@inria.fr

Alexander Lorz
email: alexan-der.lorz.1@kaust.edu.sa

Emmanuel Trélat
email: emmanuel.trelat@upmc.fr

Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy

published or not. The documents may come

Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy

Introduction

One of the primary causes of death worldwide is cancer [START_REF] Siegel | Cancer Statistics[END_REF]. Cancer treatment encounters two main pitfalls: the emergence of drug resistance in cancer cells and toxic side effects to healthy cells. Given these causes of treatment failure, designing optimized therapeutic strategies is a major objective for oncologists. In this paper, we propose a mathematical framework for modelling these phenomena and optimally combining therapies.

Overview and motivation

The most frequently used class of anti-cancer drugs are chemotherapeutic (cytotoxic) drugs, which are toxic to cells, leading to cell death. For example, platinum-based agents kill dividing cells by causing DNA damage and disrupting DNA replication [START_REF] Kelland | The resurgence of platinum-based cancer chemotherapy[END_REF]. Another class of drugs are cytostatic drugs, which slow down cell proliferation without killing cells. For example, trastuzumab is a cytostatic drug used in breast cancer treatment that targets growth factor receptors present on the surface of cells, and inhibits their proliferation [START_REF] Hortobagyi | Trastuzumab in the treatment of breast cancer[END_REF]. Despite this obvious functional difference between the two classes, cytostatic drugs, such as tyrosine kinase inhibitors, can also be cytotoxic at high doses [START_REF] Rixe | Is cell death a critical end point for anticancer therapies or is cytostasis sufficient[END_REF].

It is a well documented fact that cytotoxic agents can fail to control cancer growth and relapse [START_REF] Hanahan | Less is more, regularly: metronomic dosing of cytotoxic drugs can target tumor angiogenesis in mice[END_REF][START_REF] Pasquier | Metronomic chemotherapy: new rationale for new directions[END_REF][START_REF] Scharovsky | Metronomic chemotherapy: changing the paradigm that more is better[END_REF]. First, eradication of the tumour cell population is compromised by the emergence of drug resistance, due to intrinsic or acquired genotypic and phenotypic heterogeneity in the cancer cell population [START_REF] Bedard | Tumour heterogeneity in the clinic[END_REF][START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF][START_REF] Greaves | Evolutionary determinants of cancer[END_REF][START_REF] Navin | Inferring tumor progression from genomic heterogeneity[END_REF], because a subpopulation of resistant cells survives and proliferates, even in the presence of further treatment with identical [START_REF] Sharma | A Chromatin-Mediated Reversible Drug-Tolerant State in Cancer Cell Subpopulations[END_REF], or higher doses [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'What does not kill me strengthens me[END_REF]. Second, chemotherapeutic treatments have unwanted side effects on healthy cells, which precludes unconstrained treatment use for fear of unwanted toxicities to major organs. It is therefore a challenge for oncologists to optimally and safely treat patients with chemotherapy.

The medical objective of killing cancer cells together with preserving healthy cells from excessive toxicity is routinely translated in mathematical terms as finding the best therapeutic strategies (i.e., below some maximum tolerated dose, referred to as MTD) in order to minimise an appropriately chosen cost function. There are many works in mathematical oncology focusing on the optimal modulation of chemotherapeutic doses and schedules designed to control cancer growth, e.g. [START_REF] Agur | Optimizing chemotherapy scheduling using local search heuristics[END_REF][START_REF] Costa | Optimal chemical control of populations developing drug resistance[END_REF][START_REF] Costa | Optimal chemotherapy: a case study with drug resistance, saturation effect, and toxicity[END_REF][START_REF] Kimmel | Control theory approach to cancer chemotherapy: Benefiting from phase dependence and overcoming drug resistance[END_REF][START_REF] Ledzewicz | Analysis of models for evolving drug resistance in cancer chemotherapy[END_REF][START_REF] Ledzewicz | Optimal controls for a model with pharmacokinetics maximizing bone marrow in cancer chemotherapy[END_REF][START_REF] Ledzewicz | Optimal and suboptimal protocols for a class of mathematical models of tumor anti-angiogenesis[END_REF][START_REF] Ledzewicz | Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy[END_REF][START_REF] Swan | Optimal control analysis in the chemotherapy of igg multiple myeloma[END_REF][START_REF] Swan | Applications of optimal control theory in biomedicine[END_REF][START_REF] Swan | Role of optimal control theory in cancer chemotherapy[END_REF].

Since using ordinary differential equations (ODEs) is a common technique for modelling the temporal dynamics of cell populations, the mathematical field of optimal control applied to ODEs has emerged as an important tool to tackle such questions (see for instance [START_REF] Schättler | Optimal Control for Mathematical Models of Cancer Therapies[END_REF] for a complete presentation). In these ODE models, toxicity can either be incorporated in the cost functional as in [START_REF] Costa | Optimal chemotherapy: a case study with drug resistance, saturation effect, and toxicity[END_REF], or by adding the dynamics of the healthy cells [START_REF] Billy | Optimisation of cancer drug treatments using cell population dynamics[END_REF]. One simple, but rather coarse, paradigm used to represent drug resistance in such ODE models is by distinguishing between sensitive and resistant cancer cell subpopulations [START_REF] Costa | Optimal chemotherapy: a case study with drug resistance, saturation effect, and toxicity[END_REF][START_REF] Ledzewicz | Analysis of models for evolving drug resistance in cancer chemotherapy[END_REF]. Herein, the main tools available to obtain rigorous results are the Pontryagin maximum principle (PMP) and geometric optimal control techniques [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF][START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF][START_REF] Schättler | Geometric optimal control: theory, methods and examples[END_REF][START_REF] Trélat | Contrôle optimal: théorie & applications[END_REF].

Another paradigm used in the mathematical modelling of drug resistance relies on the idea that phenotypic heterogeneity in cancer cells and the dynamics of cancer cell populations can be understood through the principles of Darwinian evolution [START_REF] Greaves | Cancer stem cells: back to Darwin?[END_REF][START_REF] Greaves | Clonal evolution in cancer[END_REF]. Given a particular tumor micro-and macro-environment (e.g., access to oxygen, nutrients, growth factors, drug exposure), the fittest cells are selected. In the case of resistance, resistant cell subpopulations are assumed to emerge and be selected for their high levels of fitness in the presence of chemotherapeutic agents. Whether they already exist in the cell population, surviving and remaining dormant at clinically undetectable small numbers, and emerging only by natural selection, or they do not exist at all initially, but emerge as a result of an evolutionary trade-off between proliferation and development of costly survival mechanisms [START_REF] Aktipis | Life history trade-offs in cancer evolution[END_REF], likely of epigenetic nature, is still difficult to decide. The two scenarios have been studied in a modelling framework in [START_REF] Chisholm | Emergence of Drug Tolerance in Cancer Cell Populations: An Evolutionary Outcome of Selection, Nongenetic Instability, and Stress-Induced Adaptation[END_REF].

Adaptive dynamics is a branch of mathematical biology that aims at modelling Darwinian selection [START_REF] Diekmann | A beginner's guide to adaptive dynamics[END_REF][START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF]. It is thus a natural theoretical framework for the representation of phenotypic evolution in proliferating cell populations exposed to anti-cancer drugs and tumor micro-environmental factors. Non-Darwinian evolutionary principles have also been proposed to take into account drug resistance phenomena [START_REF] Pisco | Non-Darwinian dynamics in therapy-induced cancer drug resistance[END_REF]. Adaptive dynamics is amenable to modelling these principles as well. To this end, stochastic or game-theoretic points of view (see [START_REF] Champagnat | From individual stochastic processes to macroscopic models in adaptive evolution[END_REF][START_REF] Hofbauer | Adaptive dynamics and evolutionary stability[END_REF]) are standard in adaptive dynamics. Apart from ODEs, partial differential equations (PDEs) and integro-differential equations (IDEs) represent other deterministic approaches. The latter ones represent the focus of our paper. For an introduction to PDE and IDE models in adaptive dynamics, we refer the interested reader to [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF].

A common feature of these modelling techniques is that the population is structured by a trait, referred to as phenotype. The resistance level of a cell to a drug therapy is an example of such trait. Often, this variable is assumed to be continuous since it is correlated with biological characteristics, e.g., the intracellular concentration of a detoxication molecule (such as reduced glutathione), the activity of detoxifying enzymes in metabolising the administered drug, or drug efflux transporters eliminating the drug. Another possible continuous structuring variable is the ability of some cancer cells to quickly change their phenotypes (otherwise said, their intrinsic plasticity) by regulating the level of DNA methylation and/or of activity of DNA methyltransferases [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF][START_REF] Sandoval | Cancer epigenomics: beyond genomics[END_REF]. This ability is also correlated with the degree of resistance to a given drug.

To this end, a relevant modelling alternative to the binary sensitive versus resistant ODE framework (as already proposed long ago in e.g., [START_REF] Costa | Optimal chemical control of populations developing drug resistance[END_REF][START_REF] Costa | Optimal chemotherapy: a case study with drug resistance, saturation effect, and toxicity[END_REF]) consists of studying the cells at the population level using structured population dynamics. Specifically, let us denote the density of cells at time t and with phenotype x by n(t, x), with x ∈ [0, 1]. The continuous phenotype x represents an abstract level of resistance (which may be molecularly related to the activity level of an ABC transporter, or to a mean level of methylation of the DNA) to a cytotoxic drug in a cell population. Such models allow for the analysis of the asymptotic behaviour in terms of an asymptotically selected phenotype and of the total population ρ(t) := 1 0 n(t, x) dx. In the classical non-local logistic model, written as the IDE ∂n ∂t (t, x) = r(x) -d(x)ρ(t) n(t, x),

where ρ(t) := 1 0 n(t, x) dx, cells proliferate at rate r(x) and die at rate d(x)ρ(t) (the more individuals, the more competition and thus death). Such equations have well known asymptotic properties, such as convergence and concentration [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Desvillettes | On selection dynamics for continuous structured populations[END_REF][START_REF] Jabin | On selection dynamics for competitive interactions[END_REF]. For large times, ρ converges to the smallest value ρ ∞ such that r -dρ ∞ 0 on [0, 1] and n(t, •) concentrates on the set of points such that r(x) -d(x)ρ ∞ = 0. The limit is thus typically expected to be a sum of Dirac masses. This phenomenon can be interpreted as the selection of dominant traits by the environment. A common strategy for proving this asymptotic behaviour consists in showing that ρ has a bounded variation (BV ) on [0, +∞), as in [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF][START_REF] Lorz | Long-term behaviour of phenotypically structured models[END_REF].

To the best of our knowledge, general results of convergence and concentration are still elusive for systems of IDEs: the methods used in the BV framework do not seem to generalise, even in the setting that is of special interest to us, namely in the case of two competitively interacting populations of (healthy and cancer) cells. This leads to the asymptotic analysis of systems of the form

∂n 1 ∂t (t, x) = r 1 (x) -d 1 (x)I 1 (t) n 1 (t, x), ∂n 2 ∂t (t, x) = r 2 (x) -d 2 (x)I 2 (t) n 2 (t, x). (1)
The competitive coupling comes from

I 1 = a 11 ρ 1 + a 12 ρ 2 , I 2 = a 22 ρ 2 + a 21 ρ 1 with ρ i (t) = 1 0 n i (t, x) dx, i = 1, 2.
In particular, it is not clear a priori whether such interactions may or may not lead to oscillatory behaviours at the level of ρ 1 , ρ 2 . We mention the work [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] where convergence and concentration are completely characterized for a 2x2 system, where a triangular coupling structure is considered. Due to an appropriately designed Lyapunov function, the results of our paper imply that convergence and concentration hold for the model [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF].

Our goal here is also to include control terms in order to model the effect of the drugs on the proliferation and death rates. The equations we will use henceforth throughout this paper are

∂n 1 ∂t (t, x) = r 1 (x) 1 + α 1 v(t) -d 1 (x)I 1 (t) -µ 1 (x)u(t) n 1 (t, x), ∂n 2 ∂t (t, x) = r 2 (x) 1 + α 2 v(t) -d 2 (x)I 2 (t) -µ 2 (x)u(t) n 2 (t, x). (2)
Here the asymptotic analysis is more complex because of the controls u and v. On a fixed time-frame (0, T ), we will search among controls u, v in BV (0, T ), since it would not biologically feasible to impose fast varying drug infusion rates to patients.

The main difference between the approach developed in this paper and the ODE ones is that the model considered throughout is an IDE approach studied from an optimal control perspective. This model is inspired by [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF][START_REF] Lorz | Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors[END_REF] where an IDE model (possibly with an additional space variable) has been used in order to model the effect of constant doses of cytotoxic and cytostatic drug chemotherapies.

However, a mathematical proof of the failure of chemotherapeutic treatments at MTD levels to eradicate cancer cell population has not yet been formulated in the IDE modelling framework, a result that we will obtain in the present study. In addition, previous works have only considered a priori prescribed drug treatment schedules.

In this paper, we show that our model ( 2) is consistent with clinical observations on the effect of constant infusion of high doses, and we address the optimal control problem of such IDE models. Our study has a potential impact for oncologists and mathematical biologists, since it provides an accurate and robust understanding of possible optimal strategies. Up to our knowledge, this is the first time that a mathematical model and its optimal control reflect the emerging fact, observed and acknowledged recently by many clinicians, that giving maximal tolerated drug doses, even periodically, may finally be detrimental and in any case is far from being an optimal strategy in view of curing cancer.

Modelling and overview of the main results

We consider both the healthy and cancer cell populations, modelled by their respective densities n H (t, x) and n C (t, x), where the variable x ∈ [0, 1], called phenotype, represents drug resistance levels: a cell of phenotype x is highly sensitive if x is close to 0, and is highly resistant if x is close to 1. Chemotherapy is modelled by two functions of time u 1 and u 2 , representing cytotoxic and cytostatic drug infusion flows, respectively, in the two cell populations. These functions are the controls, assumed to be bounded variation functions of time if the equation is set on a bounded time-frame, subject to maximum tolerated doses (MTD) thresholds:

0 u 1 (t) u max 1 , 0 u 2 (t) u max 2 . (3)
We assume that the densities n H and n C satisfy the following lDE system:

∂n H ∂t (t, x) = R H (x, ρ H (t), ρ C (t), u 1 (t), u 2 (t)) n H (t, x), ∂n C ∂t (t, x) = R C (x, ρ C (t), ρ H (t), u 1 (t), u 2 (t)) n C (t, x), (4)
with the net growth rates defined as

R H (x, ρ H , ρ C , u 1 , u 2 ) := r H (x) 1 + α H u 2 -d H (x)I H -u 1 µ H (x), R C (x, ρ C , ρ H , u 1 , u 2 ) := r C (x) 1 + α C u 2 -d C (x)I C -u 1 µ C (x), (5)
the non-local coupling as

I H := a HH ρ H + a HC ρ C , I C := a CH ρ H + a CC ρ C , (6)
with

ρ H (t) = 1 0 n H (t, x) dx, ρ C (t) = 1 0 n C (t, x) dx,
which are the total number of healthy and tumour cells at time t.

The system starts from the initial conditions

n H (0, x) = n 0 H (x), n C (0, x) = n 0 C (x). (7)
In the above model:

• r H and r C are the drug-free proliferation rates, assumed to be positive, decreasing functions on [0, 1].

• The factors t) model the decrease in proliferation rates due to cytostatic drugs. The positive constants α H and α C represent average sensitivities of cells to cytostatic drugs. Throughout, we make the assumption that cancer cells are more sensitive to the drugs, i.e.,

1 1+α H u 2 (t) and 1 1+α C u 2 (
α H < α C . (8)
• The terms d H I H and d C I C are the drug-free death rates. The functions d H and d C are positive, decreasing functions on [0, 1]. Given the dependence in I H and I C , the model resembles a logistic one. According to their definitions, these functions are linear combinations of the total population ρ H and ρ C , i.e., we consider both intraspecific and interspecific competition. We assume that the intraspecific competition is stronger than the interspecific one:

0 < a HC < a HH , 0 < a CH < a CC . (9)
• The terms u 1 (t)µ H (x) and u 1 (t)µ C (x) are additional death rates due to cytotoxic drugs, with µ H and µ C assumed to be non-negative, decreasing functions on [0, 1]. These functions may vanish on some interval [1 -, 1], which in this case reflects the fact that some cells become fully resistant to those drugs.

• This model imposes that the phenotype x be mostly linked to the resistance to cytotoxic drugs, in accordance with [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'What does not kill me strengthens me[END_REF]. This fact which will be made clearer by the analysis of the model (see below among the consequences of Theorem 1).

Asymptotic behaviour for controls in BV ([0, +∞)). Our first aim is to show that our model reproduces the following clinical observations: when high drug doses are administered, the tumour first reduces in size before regrowing, insensitive to further treatment.

We thus want to establish asymptotic properties of the model, a challenging task since it is similar to system (2). The following statement is our first main result: we achieve a complete description of the asymptotic behaviour of system (4), with a class of asymptotically constant controls.

Theorem 1. Let u 1 , u 2 be any functions in BV ([0, +∞)), and let ū1 , ū2 be their limits at +∞. Then, for any positive initial population of healthy and of tumour cells, (ρ

H (t), ρ C (t)) converges to some equilibrium point (ρ ∞ H , ρ ∞ C
), which can be explicitly computed. Furthermore, n H and n C concentrate on a set of points which can also be explicitly computed.

The explicit values can be found in Section 2, where this result is proved. If ū1 = 0, the sets of points on which n H and n C concentrate are independent of ū2 . This is due to the fact that the phenotypic variable x models resistance to cytotoxic drugs. Thus, system (4) also exhibits convergence and concentration: the classical features of a single IDE generalise to a system, and our method is flexible enough to incorporate controls in BV ([0, +∞)). The main ingredient of the proof is an appropriately designed Lyapunov function, inspired by [START_REF] Jabin | On selection dynamics for competitive interactions[END_REF] and by the classical Lyapunov functions used for studying the global asymptotical stability of steady-states in Lotka-Volterra ODE systems, as in [START_REF] Goh | Global stability in many-species systems[END_REF]. The proof then consists of a fine analysis of this Lyapunov functional, yielding both convergence and concentration, while providing estimates on their speed. This method is new in the analysis of such IDE systems.

If µ C vanishes identically on some interval [1-ε, 1] (meaning that full resistance is possible), this theorem explains why, in the long run, high doses are not optimal. This means that our mathematical conclusions are in agreement with the idea that the standard method used in the clinic, namely administering maximum tolerated doses, should be reconsidered. Alternatives are currently extensively being investigated by oncologists, e.g., metronomic scheduling, which relies on frequent and continuous low doses of chemotherapy [START_REF] Benzekry | Maximum tolerated dose versus metronomic scheduling in the treatment of metastatic cancers[END_REF][START_REF] Carrère | Optimization of an in vitro chemotherapy to avoid resistant tumours[END_REF][START_REF] Pasquier | Metronomic chemotherapy: new rationale for new directions[END_REF].

Theorem 1 thus motivates the optimal control problem of searching for the best possible functions u 1 and u 2 to minimise the number of cancer cells within a given horizon of time, which we now introduce in more details.

Optimal control problem: optimal chemotherapy strategy. We fix some T > 0, and assume that the initial conditions n 0 H and n 0 C are continuous and positive functions on [0, 1]. For any (u 1 , u 2 ) in (BV (0, T )) 2 which satisfy (3), we consider the associated trajectory (n H (•, x), n C (•, x)) on [0, T ], solution of the system (4) starting from (n 0 H , n 0 C ). We also take into account two state constraints:

• it is required to keep a minimal proportion of healthy cells with respect to the total number of cells, and hence we impose that

ρ H (t) ρ H (t) + ρ C (t) θ HC , (10)
for some 0 < θ HC < 1.

• Moreover, we impose that the number of healthy cells always remains above a certain fraction of the initial number of healthy cells:

ρ H (t) ≥ θ H ρ H (0), (11)
for some 0 < θ H < 1.

We define A T as the set of admissible controls, i.e., for which those constraints are satisfied on (0, T ). For given (u 1 , u 2 ) ∈ A T , we define the associated cost as the number of cancer cells at the end of the time-frame:

C T (u 1 , u 2 ) := ρ C (T ). (12)
We now define the optimal control problem, denoted in short (OCP) in the sequel, as

inf (u 1 ,u 2 )∈A T C T (u 1 , u 2 ). (13)
In other words, we want to find the best drug administration strategy to minimise the number of cancer cells at the end of a fixed time-frame [0, T ], while both keeping toxicity to a tolerable level and controlling tumor size. It might seem more natural to study the problem in free final time, but as explained later on, the mapping T -→ ρ C (T ) (where ρ C (T ) is the optimal value obtained by solving (OCP) on [0, T ]) is decreasing in T . This implies that the optimal control problem in free final time T is ill-posed and does not admit any solution. The other implication is that when solving the optimal control problem in free final time t f under the constraint t f T (where T is a horizon), then the optimal solution will be such that t f = T . This is why we focus on an optimal control problem in fixed final time.

In this paper, we perform a thorough study of (OCP), both theoretically and numerically.

The theoretical study is made on a smaller class of controls which, after a long phase of constant doses, are allowed to vary on a small final time frame. More precisely, for a given T 1 < T , we consider the subclass B T ⊂ A T defined by

B T := (u 1 , u 2 ) ∈ A T , (u 1 (t), u 2 (t)) = (ū 1 , ū2 ) on (0, T 1 ), T -T 1 T M 2
where T is large and where the optimal length of the second phase T 2 := T -T 1 is bounded above by some small constant T M 2 . Thus, the first phase is long. Optimising within this class is equivalent to searching for constant optimal values ū1 , ū2 of the controls during the first phase, the length of the second phase T 2 T M 2 , and optimal BV (T 1 , T ) controls u 1 , u 2 on (T 1 , T ). The reason for this restriction to this class of controls comes from the answer to the following question: given a specific tumour size (i.e., a given number of cancer cells), what would be the optimal phenotypic cellular distribution in order tumor burden at the end of the time interval? Proposition 1 shows that, for a very short time, it is always better that the cancer cell population be concentrated on some appropriate phenotype, i.e., that the initial population be a Dirac mass at some appropriate point.

From Theorem 1, we know that it is possible to asymptotically reach Dirac masses with constant controls. The combination of these two results justifies the analysis in B T .

In this class of controls, our second main result characterises a quasi-optimal strategy in large time, a result which we now state informally.

Consider the problem of minimising ρ C (T ) within the class B T . When T is large enough, the optimal strategy approximately consists of:

• a first long-time arc, with constant controls on [0, T 1 ], at the end of which populations have almost concentrated in phenotype (for T large);

• a last short-time part, on [T 1 , T ] consisting of at most three arcs (for small T 2 = T -T 1 ):

-a boundary arc1 -a free arc with controls u 1 = u max 1 and u 2 = u max 2 ;

-a boundary arc along the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] with

u 2 = u max 2 .
The precise result and hypotheses are given by Theorem 2 in Section 3.

In order for Theorem 2 to hold, an important assumption we make is that when cancer cells are concentrated on a sensitive phenotype, the maximum tolerated doses will kill more cancer cells than healthy ones. Without this assumption, it is not clear whether one can expect the same strategy to be optimal, nor whether the patient can efficiently be treated.

We also emphasise that, for these IDEs, a PMP can be established but would not lead to tractable equations. The key property to still be able to identify the optimal strategy in B T is that the long first phase allows us to use Theorem 1: both populations concentrate and their dynamics on the last phase are (approximately) governed by ODEs, as proved in Lemma 5. The second phase can thus be analysed with ODE techniques, here the Pontryagin maximum principle (see [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF][START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]). This is done in Proposition 2.

More concretely, Theorem 2 says that:

To optimally treat a cancer, the quasi-optimal strategy consists of:

• first, administering constant doses to the patient, over a long time. The role of the first long-time arc is to allow the cancer cell population to concentrate on a sensitive phenotype. From a mathematical point of view, this means that the healthy and tumour cell populations (almost) converge to a Dirac mass.

• second, during a short-time phase, following a strategy composed of at most three arcs. If the first phase is such that the constraint (10) is saturated, then there can be a first arc along this constraint. The maximal amount of drugs is administered until the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] saturates. The last arc is along this constraint, with an appropriately chosen cytotoxic drug infusion which leads to a further decrease of the number of cancer cells.

Numerically, we solve the problem (OCP) in A T . The simulations confirm the theoretical results and show that, with the chosen set of parameters, the strategy indeed approximately consists of these two phases for T large. We also compare the optimal strategy with a periodic one, and verify that the former performs better than the latter.

Furthermore, the numerical results suggest that for generic parameters, the optimal choice of constant controls on the first phase is such that the constraint (10) is saturated. Thus, the second phase possibly starts on this constraint.

Another important property highlighted by the numerical simulations is that, given the choice of parameters made, ρ C can decrease arbitrarily close to 0 once the cancer cell population has concentrated on a sensitive enough phenotype. We thus find a strategy for which T -→ ρ C (T ) is decreasing to 0; hence, there would be no solution to (OCP) if the final time T were let free. This is the first time that a mathematical model based on integro-differential equations demonstrates that, within our modeling framework, immediate administration of maximal tolerated drug doses, or a periodic treatment schedule, is an optimal solution for eradicating cancer. Here, we prove that it is better to allow the phenotypes to concentrate, before administering maximal doses. Such a strategy is much more efficient. This is also a message to be conveyed to clinicians, who have become increasingly aware of such facts.

The paper is organised as follows. Section 2 is devoted to the proof of Theorem 1 and to numerical simulations showing how the model can reproduce the regrowth of a cancer cell population. Using these results, we have theoretical and numerical grounds for our claim that constant doses are sub-optimal and we then turn our attention to (OCP). In Section 3, several arguments are given to justify the restriction to the class B T , with a long first phase. The rest of the section is then devoted to proving Theorem 2. The numerical solutions of (OCP) in A T are provided in Section 4. They are compared to periodic strategies. In Section 5, we conclude with several comments and open questions.

Constant infusion strategies

This section is devoted to the asymptotic analysis of the IDE model ( 4), in order to specifically understand the effect of giving constant doses on the long run. We start by considering one equation only, for which relatively simple and well known arguments are sufficient to conclude that both convergence and concentration hold.

Asymptotics for healthy or cancer cells alone

In this section, we assume that n 0 H = 0 and that n 0 C is continuous and positive on [0, 1]. We have the following result of convergence and concentration for constant controls and cancer cells alone. Of course, we have a similar statement for healthy cells alone. Lemma 1. Assume that u 1 and u 2 are constant: u 1 ≡ ū1 , and u 2 ≡ ū2 , and assume that

r C 1 + α C ū2 -ū1 µ C > 0 on [0, 1]. ( 14
)
Then the total population of cancer cells ρ C (t) converges to ρ ∞ C 0, which is the smallest nonnegative real number such that

r C 1 + α C ū2 -ū1 µ C d C a CC ρ ∞ C on [0, 1]. ( 15
)
Let B C ⊂ [0, 1] be the set of points such that the equality holds in [START_REF] Chisholm | Emergence of Drug Tolerance in Cancer Cell Populations: An Evolutionary Outcome of Selection, Nongenetic Instability, and Stress-Induced Adaptation[END_REF].

Then n C (t, •) concentrates on B C as t goes to +∞. In particular, if B C is reduced to a singleton x ∞ C , then n C (t, •) converges to ρ ∞ C δ x ∞ C in M 1 (0, 1).
Here and in the sequel, δ x denotes the Dirac mass at x, and M 1 (0, 1) is the set of Radon measures supported in [0,1].

The proof of this lemma is rather classical, its main ingredient is proving that ρ is a BV function to obtain convergence, and concentration follows. It is done in Appendix A. Remark 1. As it clearly appears in the proof, the result holds for more general initial conditions in L ∞ (0, 1). One only needs to require that they are bounded from below by a positive constant on a neighbourhood of one of the points of B C .

Asymptotics for the complete model: proof of Theorem 1

Now, let us take into account the complete coupling between healthy and tumour cells. For the remaining part of this article, we assume for simplicity that both n 0 H and n 0 C are continuous and positive on [0, 1] (although it is possible to be slightly more general, see Remark 1), but we emphasise that we no longer require assumption [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF]. A further technical assumption is needed to prove that convergence and concentration hold, namely that the functions are Lipschitz continuous:

r H , r C , d H , d C , µ H , µ C ∈ C 0,1 (0, 1). (16)
In two dimensions and with constant controls ū1 , ū2 , the previous technique of proving that ρ is BV cannot be extended. As for a single equation, however, we can integrate the equations with respect to x to obtain upper bounds for ρ H and ρ C . For example, let us integrate the equation defining ρ H and bound as follows:

dρ H (t) dt 1 0 r H (x) -d H (x)a HH ρ H (t) n H (t, x) dx.
Thus, we clearly have lim sup t→+∞ ρ i ρ max i := max x r i a ii d i for i = H, C. It also still holds with the reasoning made in the proof of Lemma 1 that if ρ H and ρ C converge, then the limits must be the solution of the (invertible) system

a HH ρ ∞ H + a HC ρ ∞ C = I ∞ H , a CH ρ ∞ H + a CC ρ ∞ C = I ∞ C . (17)
where

I ∞ H 0 is the smallest nonnegative real number such that r H (x) 1 + α H ū2 -ū1 µ H (x) d H (x)I ∞ H , (18)
and I ∞ C 0 is the smallest nonnegative real number such that

r C (x) 1 + α C ū2 -ū1 µ C (x) d C (x)I ∞ C . (19)
Furthermore, if this convergence holds true, then n H (resp. n C ) concentrate on A H (resp. A C ) defined as

A H = x ∈ [0, 1], r H (x) 1 + α H ū2 -ū1 µ H (x) -d H (x)I ∞ H = 0 , A C = x ∈ [0, 1], r C (x) 1 + α C ū2 -ū1 µ C (x) -d C (x)I ∞ C = 0 .
Proof of Theorem 1.

First step: definition of the Lyapunov functional.

We adapt a strategy developed in [START_REF] Jabin | On selection dynamics for competitive interactions[END_REF]. We choose any couple of measures

(n ∞ H , n ∞ C ) in M 1 (0, 1) satisfying 1 0 n ∞ i (x) dx = ρ ∞ i , i = H, C which furthermore satisfy supp(n ∞ H ) ⊂ A H , supp(n ∞ C ) ⊂ A C . (20)
For i = H, C, and m i := 1 d i , let us define the Lyapunov functional as

V (t) := λ H V H (t) + λ C V C (t),
where

V i (t) = 1 0 m i (x) n ∞ i (x) ln 1 n i (t, x) + (n i (t, x) -n ∞ i (x)) dx,
with positive constants λ H and λ C to be adequately chosen later.

Second step: computation and sign of the derivative.

In what follows, we skip dependence in t in the functions R H and R C to increase readability. We have

dV H dt = 1 0 m H (x) -n ∞ H (x) ∂ t n H (t, x) n H (t, x) + ∂ t n H (t, x) dx = 1 0 m H (x) R H (x, ρ H , ρ C , u 1 , u 2 ) [n H (t, x) -n ∞ H (x)] dx = 1 0 m H (x) (R H (x, ρ H , ρ C , u 1 , u 2 ) -R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 )) [n H (t, x) -n ∞ H (x)] dx + 1 0 m H (x) R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 ) [n H (t, x) -n ∞ H (x)] dx
The first term is simply

1 0 m H (x) (R H (x, ρ H , ρ C , u 1 , u 2 ) -R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 )) [n H (t, x) -n ∞ H (x)] = 1 0 m H (x)d H (x) [a HH (ρ ∞ H -ρ H ) + a HC (ρ ∞ C -ρ C )] [n H (t, x) -n ∞ H (x)] dx = -a HH (ρ ∞ H -ρ H ) 2 -a HC (ρ ∞ C -ρ C )(ρ ∞ H -ρ H )
The second term can also be written as

B H (t) := 1 0 m H (x) R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 ) [n H (t, x) -n ∞ H (x)] dx = 1 0 m H (x) R H (x, ρ ∞ H , ρ ∞ C , ū1 , ū2 ) [n H (t, x) -n ∞ H (x)] dx + 1 0 m H (x) (R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 ) -R H (x, ρ ∞ H , ρ ∞ C , ū1 , ū2 )) [n H (t, x) -n ∞ H (x)] dx = 1 0 m H (x)R H (x, ρ ∞ H , ρ ∞ C , ū1 , ū2 ) n H (t, x) dx + 1 0 m H (x) r H (x) 1 1 + α H u 2 - 1 1 + α H ū2 + µ H (x)(ū 1 -u 1 ) [n H (t, x) -n ∞ H (x)] dx,
where we use [START_REF] Desvillettes | On selection dynamics for continuous structured populations[END_REF] for the last equality. Note that the first term in the last expression is nonpositive by definition of (ρ ∞ H , ρ ∞ C ), and the second goes to 0 as t goes to +∞. Consequently, the decomposition

B i = Bi + E i , i = H, C, (21)
holds, with BH , BC nonpositive, and E H , E C which asymptotically vanish. This decomposition will be important in the last step.

Eventually, we have:

dV dt = - 1 2 X T M X + λ H B H + λ C B C (22)
with

M = A T D + DA, X = ρ ∞ H -ρ H ρ ∞ C -ρ C , D = λ H 0 0 λ C and A = a HH a HC a CH a CC .
We first look for a choice of constants λ H , λ C that ensures that the symmetric matrix M , is also positive semi-definite. using the assumption (9).

Since M = 2λ H a HH λ H a HC + λ C a CH λ H a HC + λ C a CH 2λ C a CC has
Our aim is to prove that -1 2 X T M X converges to 0 as t goes to +∞, which will yield the convergence of (ρ H , ρ C ). Concentration of (n H , n C ) then follows easily with the arguments developed in the proof of Lemma 1 .

Third step: lower estimate for V.

To estimate V from below, we need a uniform (in x) upper bound on n H , n C . Because of the regularity assumption made on the data (functions are Lipschitz continuous), there exists C > 0 such that:

∀(x, y) ∈ [0, 1], R H (y, ρ H , ρ C , u 1 , u 2 ) R H (x, ρ H , ρ C , u 1 , u 2 ) -C|x -y|,
C can be chosen to be independent of t since ρ H , ρ C , u 1 , and u 2 are all bounded. This implies that

1 0 n H (t, y) dy = 1 0 n 0 H (y) exp t 0 R H (y, ρ H , ρ C , u 1 , u 2 ) ds dy 1 0 n 0 H (y) n 0 H (x) n 0 H (x) exp t 0 R H (x, ρ H , ρ C , u 1 , u 2 ) exp (-Ct|x -y|) dy n H (t, x) n 0 H (x) 1 0 exp (-Ct|x -y|) dy
Using the boundedness of ρ H and n 0 H (C has changed and is independent of t and x) and computing the integral, we can write: for t large enough, n H (t, x)

Ct. Similarly, n C (t, x) Ct. The bound on V follows immediately:

V (t) -C (ln(t) + 1) . (23)
We now define another function, close to dV dt , whose behaviour will allow us to conclude.

Fourth step: estimates on dV dt . We set

G := - 1 2 X T M X + 2 (λ H B H + λ C B C ) .
The first term is:

- 1 2 X T M X = -λ H a HH (ρ ∞ H -ρ H ) 2 -λ C a CC (ρ ∞ C -ρ C ) 2 -(λ H a HC + λ C a CH )(ρ ∞ H -ρ H )(ρ ∞ C -ρ C ), so that (writing in short R H for R H (x, ρ H , ρ C , u 1 , u 2 ) and R C for R C (x, ρ C , ρ H , u 1 , u 2 )): d dt - 1 2 X T M X = (ρ ∞ H -ρ H ) 2λ H a HH 1 0 R H n H + (λ H a HC + λ C a CH ) 1 0 R C n C + (ρ ∞ C -ρ C ) 2λ C a CC 1 0 R C n C + (λ H a HC + λ C a CH ) 1 0 R H n H
Let us now differentiate B H and B C using the expression that initially defines them:

dB H dt = 1 0 m H (x)R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 )R H (x, ρ H , ρ C , u 1 , u 2 )n H (t, x) dx - 1 0 m H (x) α H r H (x) 1 (1 + α H u 2 ) 2 du 2 dt + µ H (x) du 1 dt [n H (t, x) -n ∞ H (x)] dx
Using the following estimate

m H (x)R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 ) R H (x, ρ H , ρ C , u 1 , u 2 ) m H (x) (R H (x, ρ ∞ H , ρ ∞ C , u 1 , u 2 ) -R H (x, ρ H , ρ C , u 1 , u 2 )) R H (x, ρ H , ρ C , u 1 , u 2 ) = -[a HH (ρ ∞ H -ρ H ) + a HC (ρ ∞ C -ρ C )] R H (x, ρ H , ρ C , u 1 , u 2 ),
we obtain the inequality

dB H dt -λ H [a HH (ρ ∞ H -ρ H ) + a HC (ρ ∞ C -ρ C )] 1 0 R H n H - 1 0 m H (x) r H (x) 1 (1 + α H u 2 ) 2 du 2 dt + µ H (x) du 1 dt [n H (t, x) -n ∞ H (x)] dx
and similarly for B C . We thus obtain (the terms in a HH and a CC cancel each other out):

dG dt (ρ ∞ H -ρ H ) (λ H a HC + λ C a CH -2λ C a CH ) 1 0 R C n C + (ρ ∞ C -ρ C ) (λ H a HC + λ C a CH -2λ H a HC ) 1 0 R H n H -a(t) du 2 dt + b(t) du 1 dt (24)
where a and b are bounded functions defined by

a(t) := 2λ H 1 (1 + α H u 2 ) 2 1 0 m H (x)r H (x)(n H (t, x) -n ∞ H (x)) dx + 2λ C 1 (1 + α C u 2 ) 2 1 0 m C (x)r C (x)(n C (t, x) -n ∞ C (x)) dx, b(t) := 2λ H 1 0 m H (x)µ H (x)(n H (t, x) -n ∞ H (x)) dx + 2λ C 1 0 m C (x)µ C (x)(n C (t, x) -n ∞ C (x)) dx.
The first two terms at the right-hand side of inequality [START_REF] Fourer | A modeling language for mathematical programming[END_REF] are equal to 0 thanks to the choice made for (λ H , λ C ) so that it simplifies to

dG dt -a(t) du 2 dt + b(t) du 1 dt . ( 25
)
Fifth step: conclusion.

Noting that dV dt 1 2 G, it follows that V (t) -V (0) 1 2 t 0 G(s) ds. From G(s) = G(t) - t s dG dt (
z) dz, using [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF] and by integrating the previous inequality, we have

V (t) -V (0) t 1 2 G(t) + 1 2t t 0 t s a(z) du 2 dt (z) + b(z) du 1 dt (z) dz ds.
Now, using the decomposition [START_REF] Diekmann | A beginner's guide to adaptive dynamics[END_REF] introduced in the second step, we obtain:

2 V (t) -V (0) t - 1 t t 0 t s a(z) du 2 dt (z) + b(z) du 1 dt (z) dz ds -2 (λ H E H + λ C E C ) - 1 2 X T M X + 2 λ H BH + λ C BC .
In other words, since the right-hand side of this inequality consists of nonpositive terms, the claim on the convergence of ρ H and ρ C is proved if we establish that the left-hand side tends to 0.

As a consequence of the estimate ( 23) on V established in the third step, 2 V (t) -V (0) t converges to 0. This is also true for 2 (λ H E H + λ C E C ). It thus remains to analyse the asymptotic behaviour of the function

1 t t 0 t s a(z) du 2 dt (z) + b(z) du 1 dt (z) dz ds.
The analysis relies on the following lemma. Lemma 2. Let φ in L ∞ (0, +∞), and u in BV ([0, +∞)). Then

lim t→+∞ 1 t t 0 t s φ(z)u (z) dz ds = 0.
Proof. Let us start by writing

1 t t 0 t s φ(z)u (z) dz ds = 1 t t 0 t 0 φ(z)u (z) dz ds - 1 t t 0 s 0 φ(z)u (z) dz ds = t 0 φ(z)u (z) dz - 1 t t 0 s 0 φ(z)u (z) dz ds = Γ(t) - 1 t t 0 Γ(s) ds
where Γ(t) := t 0 φ(z)u (z) dz. The expression above can thus be decomposed as the function Γ minus its Cesàro average. To conclude, it suffices that Γ has a limit at +∞, which in turn is true as soon as φ u is integrable on the half-line. This fact is a direct consequence of the boundedness of φ and the integrability of the derivative of a BV function on [0, +∞).

This ends the proof of Theorem 1. Remark 2. The situation differs from [START_REF] Jabin | On selection dynamics for competitive interactions[END_REF] where the non-local logistic term is of the form 1 0 b(x, y) n(t, y) dy, with some strong competition assumption on the kernel b. In particular, that assumption implies the uniqueness of the ESDs, which is not necessarily true in our setting. Remark 3. This theorem means that under general conditions, both populations concentrate and the total number of healthy (resp., cancer) cells converge. In the case of constant controls and when there is selection of a unique phenotype in both populations, it provides a complete understanding of the mapping

(ū 1 , ū2 ) -→ (x ∞ H , x ∞ C , ρ ∞ H , ρ ∞ C ), (26)
where

ρ ∞ H δ x ∞ H and ρ ∞ C δ x ∞
C are the respective limits of n H (t, •) and n C (t, •) in M 1 (0, 1), as t goes to +∞. In particular, if we restrict ourselves to constant controls and a large time T , the problem of minimising ρ C (T ) is equivalent to minimising ρ ∞ C as a function of (ū 1 , ū2 ).

Speed of convergence in Theorem 1

Because it relies on a Lyapunov functional, Theorem 1 also yields results on the speeds of convergence when the controls are constant, which we state and analyse separately in the following corollary.

Corollary 1. Assume u 1 ≡ ū1 , u 2 ≡ ū2 .
Then, as t → +∞, (ρ

H , ρ C ) = (ρ ∞ H , ρ ∞ C ) + O ln(t) t 1 2
and concentration occurs at speed

O ln(t) t
, in the following sense:

1 0 m H (x)R H (x, ρ ∞ H , ρ ∞ C , ū1 , ū2 ) n H (t, x) dx = O ln(t) t , 1 0 m C (x)R C (x, ρ ∞ C , ρ ∞ H , ū1 , ū2 ) n C (t, x) dx = O ln(t) t .
In particular, if A H is reduced to a singleton x ∞ H , then

∀ > 0, [0,1]\[x ∞ H -,x ∞ H + ] n H (t, x) dx = O ln(t) t ,
and similarly for A C .

Proof. The speed of convergence of (ρ H , ρ C ) and of the integrals can be obtained by rewriting the end of the proof, namely that -1 2 X T M X + 2 λ H BH + λ C BC is bounded from below by 2 V (t)-V (0) t , which converges to 0 as O ln(t) t . Since each of those three terms is nonpositive, they all converge to 0 at the previous speed, and the integrals of interest are nothing but the functions BH and BC .

For the last statement, we fix > 0 and denote h

:= -m H R H (•, ρ ∞ H , ρ ∞ C , ū1 , ū2 ) 0 on [0, 1],
which by assumption vanishes at x ∞ H only. We choose a > 0 small enough such that

a1 [0,1]\[x ∞ H -,x ∞ H + ]
h on [0, 1]. This enables us to write

[0,1]\[x ∞ H -,x ∞ H + ] n H (t, x) dx 1 a 1 0 m H (x)R H (x, ρ ∞ H , ρ ∞ C , ū1 , ū2 ) n H (t, x) dx = O ln(t) t .
The reasoning is the same if A C is reduced to a singleton.

Remark 4. Although the Lyapunov functional gives us information on the speed of both phenomena in the sense defined above, it does not say whether one of the two is faster.

However, if A H or A C is reduced to a singleton, the speed O ln(t) t obtained for the convergence to 0 of the expression -1 2 X T M X + 2 λ H BH + λ C BC is almost optimal: there cannot exist any α > 1 such that this sum vanishes like O 1 t α . This comes from the fact that if it were to hold true, dV dt would be integrable on the half-line, which would imply the convergence of V . This is not possible since either V H or V C goes to -∞.

Mathematical simulations of the effect of constant drug doses

Throughout the study, we will consider the following numerical data, taken from [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF]:

r H (x) = 1.5 1 + x 2 , r C (x) = 3 1 + x 2 , d H (x) = 1 2 (1 -0.1x), d C (x) = 1 2 (1 -0.3x), a HH = 1, a CC = 1, a HC = 0.07, a CH = 0.01, α H = 0.01, α C = 1, u max 1 = 3.5, u max 2 = 7,
and the initial data

n H (0, x) = K H,0 exp(-(x -0.5) 2 /ε), n C (0, x) = K C,0 exp(-(x -0.5) 2 /ε),
with ε > 0 small (typically, we will take either ε = 0.1 or ε = 0.01), and where K H,0 > 0 and K C,0 > 0 are such that ρ H (0) = 2.7, ρ C (0) = 0.5.

The value ρ H (0) is not the same as in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF]: it is chosen to be slightly below the equilibrium value of the system with n C ≡ 0, u 1 ≡ 0, u 2 ≡ 0, in accordance with the fact that there is homeostasis in a healthy tissue. Indeed, we start with a non-negligible tumour which must have (due to competition) slightly lowered the number of healthy cells with comparison to a normal situation.

We also define ρ CS (t) :

= 1 0 (1 -x)n C (t,
x) dx, which may be seen as the total number at time t of tumour cells that are sensitive, and ρ CR (t) := 1 0 xn C (t, x) dx, which may be seen as the total number at time t of tumour cells that are resistant.

Of course, sensitivity/resistance being by construction a non-binary variable, the weights x and 1 -x are an example of a partition between a relatively sensitive class and a relatively resistant class in the cancer cell population; other choices might be made for these weights, e.g., x 2 and 1 -x 2 .

Discussion of the choice for µ H and µ C . These functions measure the efficiency of the drugs treatment. The choice done in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF] is

µ H (x) = 0.2 0.7 2 + x 2 , µ C (x) = 0.4 0.7 2 + x 2 .
However, with this choice of functions, if we take constant controls u 1 and u 2 , with

u 1 (t) = Cst = u max 1 = 3.5, u 2 (t) = Cst = 2,
then we can kill all tumour cells (at least, they decrease exponentially to 0), and no optimisation is necessary. The results of a simulation can be seen on Figure 1. ρ H (t)+ρ C (t) and of t → ρ CS (t) ρ C (t) .

On Figure 1, the population of tumour cells is Gaussian-shaped, decreases exponentially to 0 while its center is being shifted to the right: it means that tumour cells become more and more resistant as time goes by. This is in agreement with the fact that cells acquire resistance to treatment when drugs are given constantly. However, although the proportion of sensitive cells t → ρ CS (t) ρ C (t) is quickly decreasing, the drugs are still efficient at killing the cells. This is not realistic, as it does not match the clinically observed saturation phenomenon. Most cancer cells have acquired resistance and any immediate further treatment should have no effect. In the simulation above, there is no saturation because the function µ C is continuous and positive over the whole interval [0, 1] and is not small enough close to 1. In order to model this saturation phenomenon, we choose to modify the model used in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], by modifying slightly the function µ C . The new function µ C that will throughout be considered is defined by µ C (x) = max 0.9 0.7 2 + 0.6x 2 -1, 0 .

On Figure 2, the former function µ C is in blue, and the new one is in red. With this new function, the simulation of Figure 1, with u 1 = Cst = 3.5 and u 2 = Cst = 2, is completely modified, as can be seen on Figure 3. Indeed, this time, the strategy consisting of taking constant controls u 1 (t) = Cst = 3.5 and u 2 (t) = Cst = 2 is not efficient anymore and does not allow for (almost) total eradication of the tumour. In sharp contrast, we observe on Figure 1 that the tumour cells are growing again, moreover concentrating around some resistant phenotype. C . However, it is quite clear that choosing the optimal constant dose (ū 1 , ū2 ) to minimise ρ ∞ C leaves room for improvement, as the choice (u max 1 , u max

2

) is still the optimal one for sensitive cells. Therefore, it makes sense to allow u 1 and u 2 to be any functions satisfying (4) as in (OCP), which we will study both from the theoretical and numerical points of view in the next two sections.

Theoretical analysis of (OCP)

Before analysing (OCP), let us first consider a much simpler ODE model for which we can find the solution explicitly in order to develop some intuition.

Simplified optimal control problems

We consider the ODE

dρ dt = (r -dρ(t) -µu(t))ρ(t), ρ(0) = ρ 0 > 0. (27)
(C1) Optimal control problem: minimise ρ(T ) over all possible solutions of ( 27) with a L 1 -constraint on u, i.e.,

T 0 u(t) dt ≤ u 1,max . (28)
Lemma 3. The optimal solution for problem (C1) is

u opt = u 1,max δ t=T .
Remark 5. The statement can be misleading: we actually prove that there is no optimal solution, but rather that the problem leads to an infimum. Still, writing u opt = u 1,max δ t=T makes sense as a way to obtain the infimum is to take a family (u ) >0 in L 1 which converges to u opt , for example u :

= 1 u 1,max 1 [T -,T ] .
Adding another the constraint, we have a second optimal control problem (C2) minimise ρ(T ) over all possible solutions of ( 27) with the L 1 -constraint ( 28) and a

L ∞ -constraint u ≤ u ∞,max . (29)
We assume u ∞,max T > u 1,max , since otherwise it is clear that the optimal strategy is u ∞,max 1 [0,T ] . Lemma 4. We define T 1 (T ) := T -u 1,max u ∞,max . The optimal solution for problem (C2) is

u op = u ∞,max 1 [T 1 ,T ] . (30)
The proofs of these two results can be found in Appendix B.

Remark 6. The previous lemmas on simplified equations give some insight on two important features:

• for large times, constant controls lead to concentration, as evidenced by Theorem 1. As explained more rigorously further below in Lemma 5, when the populations are concentrated on some single phenotypes, the integro-differential equations boil down to ODEs, for which the last results and standard techniques from optimal control theory apply.

• for ODE models, it is optimal to use the maximal amount of drug at the end of the timewindow if there is a L 1 constraint on the control. Avoiding the emergence of resistance will indirectly act as some L 1 constraint, which is why this result also provides some interesting intuition on the optimal control problem (OCP).

Assumptions and further remarks

Let us start by mentioning a possible alternative state constraint for (OCP).

Remark 7. Alternatively to [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF], we might want to directly control the number of cancer cells and replace [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF] by

ρ C (t) C max (31)
for some C max > 0. The set of constraints ( 10), [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] on the one hand, and ( 31), [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] on the other hand, are similar. Although we focus on the first one in the sequel, our analysis applies to the other set of constraints.

We now make several important additional assumptions which will be used throughout this section, all relying on the notations of Theorem 1. Our first assumption is that

for any 0 ū1 u max 1 , 0 ū2 u max 2
, A H and A C are reduced to singletons.

In this case, recall that Theorem 1 provides a mapping (ū

1 , ū2 ) → (x ∞ H , x ∞ C , ρ ∞ H , ρ ∞ C )
, and with a slight abuse of notation, we will omit the dependence in (ū 1 , ū2 ) in the following final assumptions: whenever ū1 , ū2 are admissible (i.e., such that neither the constraint (10) nor the constraint ( 11) is violated), we require that the solution of the ODE system

dρ H dt = R H (x ∞ H , ρ H , ρ C , u max 1 , u max 2 ) ρ H , dρ C dt = R C (x ∞ C , ρ C , ρ H , u max 1 , u max 2 ) ρ C , (33)
with initial data (ρ ∞ C , ρ ∞ H ), has the following properties:

d dt ρ C < 0, d dt ρ H < 0 and d dt ρ C ρ H < 0. (34)
The assumption [START_REF] Jabin | On selection dynamics for competitive interactions[END_REF] means that both populations of cells decrease but that the treatment is more efficient on cancer cells. In some sense, this is a curability assumption and it will be crucial in the sequel.

We now motivate the choice of restricting our attention to the class B T by giving two results.

Optimality of a concentrated initial population for a small time

Here, we assume that for any 0

ρ C ρ max C , 0 ρ H ρ max H , 0 u 1 u max 1 and 0 u 2 u max 2 , x → R C (x, ρ C , ρ H , u 1 , u 2 ) has a unique minimum. ( 35
)
For a given initial amount of cancer cells ρ 0 C > 0, we define:

A ρ 0 C := n 0 C ∈ M 1 (0, 1) such that 1 0 n 0 C (x) dx = ρ 0 C . For n 0 C ∈ A ρ 0 C
, and given n 0 H in M 1 (0, 1), final time t f > 0, and controls u 1 , u 2 in BV (0, t f ) satisfying (3), we consider the associated trajectory (n H (•, x), n C (•, x)) on [0, t f ] solution to the system (4) starting from (n 0 H , n 0 C ). We consider the following minimisation problem

inf 0 u 1 (t) u max 1 0 u 2 (t) u max 2 inf n 0 C ∈Aρ C 0 ρ C (t f ). ( 36
)
In other words, for a fixed initial tumour size, we aim at tackling the following question:

what is the cancer cells' best possible repartition in phenotype?

A simpler (and instantaneous) version of the previous optimisation problem for t f small is inf

0 u 1 u max 1 0 u 2 u max 2 inf n 0 C ∈Aρ C 0 dρ C dt (0), (37)
for which the solution is easily obtained, and given in the following proposition.

Proposition 1. Let g := R C •, ρ 0 C , ρ 0 H , u max 1 , u max 2
. We define x C by {x C } := arg min g and ñ0

C := ρ 0 C δ x C . The optimal solution for the optimisation problem ( 37) is given by

u max 1 , u max 2 , ρ 0 C δ x C . (38)
Proof. For any

0 u 1 u max 1 , 0 u 2 u max 2 , n 0 C ∈ A ρ C 0 , dρ C dt (0) = 1 0 R C x, ρ 0 C , ρ 0 H , u 1 , u 2 n 0 C (x) dx 1 0 g(x) n 0 C (x) dx
with equality if and only if

u 1 = u max 1 , u 2 = u max 2 .
We also have

1 0 g(x) n 0 C (x) dx 1 0 g(x C ) n 0 C (x) dx = g(x C )ρ 0
C and it remains to prove that there is equality if and only if

n 0 C = ρ 0 C δ x C . If n 0 C = ρ 0 C δ x C there exists a ∈ supp n 0 C , a = x C : it is therefore possible to find > 0 such that both x C ∈ [a -, a + ] and [a-,a+ ] n 0 C (x) dx > 0. This implies 1 0 (g(x) -g(x C )) n 0 C (x) dx [a-,a+ ] (g(x) -g(x C )) n 0 C (x) dx > 0,
which concludes the proof.

For (OCP), the previous Proposition means that, very close to T , the best shape of the cancer cell density n C (t, •) is a Dirac mass. As it was proved in Theorem 1, it is possible (in arbitrarily large time) to reach Dirac masses with constant controls. The combination of these two results is our motivation for the restriction to the set B T .

Reduction of IDEs to ODEs at the end of the long first phase

Because of the previous result, it makes sense to steer the cancer dell density as close as possible to a Dirac mass. As it was proved in Theorem 1, it is possible (in large time limit) to reach Dirac masses with constant controls. Our aim is now to prove that if we give constant controls (ū 1 , ū2 ) for a long time, the dynamics of the total number of cells (ρ H , ρ C ) are arbitrarily close to being driven by a system of ODEs, a result which comes from the concentration of the IDE on (x ∞ H , x ∞ C ). The rigorous statement is given hereafter: Lemma 5. We fix T 2 > 0, 0 ū1 u max 1 and 0 ū2 u max 2 . We consider any controls (u 1 , u 2 ) defined on [0, T 1 + T 2 ] as follows: they are constant equal to (ū 1 , ū2 ) on [0, T 1 ], and any BV functions on [T 1 , T 1 + T 2 ] which satisfy [START_REF] Aktipis | Life history trade-offs in cancer evolution[END_REF]. Let (n H , n C ) be the solution of ( 4) on [0, T 1 + T 2 ], with corresponding (ρ H , ρ C ).

Then lim

T 1 →+∞ sup [T 1 ,T 1 +T 2 ] max (|ρ H -ρH |, |ρ C -ρC |) = 0,
where (ρ H , ρC ) solves the controlled ODE system

dρ H dt = R H (x ∞ H , ρH , ρC , u 1 , u 2 ) ρH , dρ C dt = R C (x ∞ C , ρC , ρH , u 1 , u 2 ) ρC , ( 39
)
defined on [T 1 , T 1 + T 2 ], starting at T 1 from (ρ H (T 1 ), ρ C (T 1 )).
Proof. Let > 0. We focus on the equation on n H which we integrate in x for any t ∈ [T 1 , T 1 + T 2 ]:

dρ H dt = 1 0 R H (x, ρ H , ρ C , u 1 , u 2 )n H (t, x) dx = R H (x ∞ H , ρ H , ρ C , u 1 , u 2 ) ρ H + 1 0 (R H (x, ρ H , ρ C , u 1 , u 2 ) -R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )) n H (t, x) dx
For the first term, we write

R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )ρ H = R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )ρ H + R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )(ρ H -ρH ) = dρ H dt + ρH d H (x ∞ H ) (-a HH (ρ H -ρH ) -a HC (ρ C -ρC )) + R H (x ∞ H , ρ H , ρ C , u 1 , u 2 ) (ρ H -ρH ) This means we end up with d dt (ρ H -ρH ) = ρH d H (x ∞ H ) (-a HH (ρ H -ρH ) -a HC (ρ C -ρC )) + R H (x ∞ H , ρ H , ρ C , u 1 , u 2 ) (ρ H -ρH ) + 1 0 (R H (x, ρ H , ρ C , u 1 , u 2 ) -R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )) n H (t, x) dx.
We look at the last term separately: the first two ones are linked to the discrepancy between ρ and ρ, while the last one will be small because n H is concentrated if T 1 is large enough. Setting w := max (|ρ H -ρH |, |ρ C -ρC |), we have the differential inequality

d dt |ρ H -ρH | Cw + 1 0 (R H (x, ρ H , ρ C , u 1 , u 2 ) -R H (x ∞ H , ρ H , ρ C , u 1 , u 2 )) n H (t, x) dx (40
) for some constant C > 0. The last term can be decomposed as

1 1 + α H u 2 1 0 (r H (x) -r H (x ∞ H )) n H (t, x) dx -u 1 1 0 (µ H (x) -µ H (x ∞ H )) n H (t, x) dx -I H 1 0 (d H (x) -d H (x ∞ H )) n H (t, x) dx.
Note that u 1 ,

1 1+α H u 2 and I H are all bounded on [T 1 , T 1 + T 2 ]
. Thus, if for any generic function φ,

1 0 (φ H (x) -φ H (x ∞ H )) n H (t,
x) dx is arbitrarily small, so is the last quantity. To that end, we write the solution of the IDE in exponential form

n H (t, x) = n H (T 1 , x) exp t T 1 R H (x, ρ H (s), ρ C (s), u 1 (s), u 2 (s)) ds ,
where the exponential is uniformly bounded on

[0, 1] × [T 1 , T 1 + T 2 ], which means that 1 0 (φ H (x) -φ H (x ∞ H )) n H (t, x) dx C 1 0 |φ H (x) -φ H (x ∞ H )| n H (T 1 , x) dx. Since n H (T 1 , •) converges to ρ ∞ H δ x ∞ H in M 1 (0, 1)
as T 1 goes to +∞, this quantity is arbitrarily small. Plugging this estimate into [START_REF] Ledzewicz | Optimal and suboptimal protocols for a class of mathematical models of tumor anti-angiogenesis[END_REF] and writing a similar inequality for the equations on the cancer cells, we obtain for T 1 large enough dw dt Cw + . We conclude by applying the Gronwall lemma, together with the fact that w(T 1 ) = 0.

Analysis of the second phase

According to the previous results, for large T and admissible constant controls (ū 1 , ū2 ), we arrive at concentrated populations whose dynamics are driven by a system of ODEs. This naturally leads to considering the following optimal control problem, on the resulting ODE concentrated in

(x ∞ H , x ∞ C ), starting from (ρ ∞ H , ρ ∞ C ) at t = 0. For readability, we write g H for g H (x ∞ H ) (resp., g C for g(x ∞ C ))
for any function g H (resp., g C ), and we stress that all assumptions made in this subsection are made for all possible admissible constant controls (ū 1 , ū2 ).

The ODE system of equations now reads

dρ H dt = r H 1 + α H u 2 -d H I H -u 1 µ H R H ρ H , dρ C dt = r C 1 + α C u 2 -d C I C -u 1 µ C R C ρ C . (41)
For a given T M 2 > 0, we investigate the optimal problem of minimising ρ C (t f ) for t f T M 2 and controls (u 1 , u 2 ) which satisfy (3), as well as the constraints [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] and [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF]. The constraint [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF]

rewrites ρ C ρ H γ with γ := 1 -θ HC θ HC .
Assume that there exists an optimal solution which is the concatenation of free and constrained arcs (either on the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] or ( 10)), with associated times (t i ) 1 i M . In particular, we thus assume without loss of generality that the parameters are such that both constraints do not saturate simultaneously on an optimal arc.

Then, by the Pontryagin maximum principle for an optimal control problem with state constraints (see [START_REF] Vinter | Optimal control. Systems and control : foundations & applications[END_REF]), there exists a bounded variation adjoint vector p = (p H , p C ) defined on [0, t f ], a scalar p 0 0, non-negative functions η 1 and η 2 and non-negative scalars ν i , i = 1, . . . , M such that if we define the Hamiltonian function by

H(ρ H , ρ C , p H , p C , u 1 , u 2 ) := p H R H ρ H + p C R C ρ C + η 1 (θ H ρ 0 H -ρ H ) + η 2 (ρ C -γρ H ) = -p H d H I H ρ H -p C d C I C ρ C + r H p H ρ H 1 + α H u 2 + r C p C ρ C 1 + α C u 2 -(µ H p H ρ H + µ C p C ρ C ) u 1 + η 1 (θ H ρ 0 H -ρ H ) + η 2 (ρ C -γρ H ),
we have 1. p, p 0 , η 1 , η 2 and the (ν i ) i=1,...,M are not all zero.

The adjoint vector satisfies

dp H dt = - ∂H ∂ρ H = -p H (-a HH d H ρ H + R H ) + a CH d C p C ρ C + η 1 + γη 2 , dp C dt = - ∂H ∂ρ C = -p C (-a CC d C ρ C + R C ) + a HC d H p H ρ H -η 2 , (43)
with

p H (t f ) = 0, p C (t f ) = p 0 . 3. t -→ η 1 (t) (resp. t -→ η 2 (t)
) is continuous along (11) (resp. ( 10)), and is such that

η 1 (θ H ρ 0 H -ρ H ) = 0 (resp. η 2 (ρ C -γρ H ) = 0) on [0, t f ]. 4.
For any i = 1, . . . , M , the Hamiltonian is continuous at t i . If t i is a junction or contact2 point with the boundary (11) (resp. with the boundary [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF]

), p H (t + i ) = p H (t - i ) + ν i , p C (t + i ) = p C (t - i ) (resp. p H (t + i ) = p H (t - i ) + γν i , p C (t + i ) = p C (t - i ) -ν i ). 5.
The controls u 1 , u 2 maximise the Hamiltonian almost everywhere.

We now make several technical assumptions (for all admissible constant controls (ū 1 , ū2 )) by requiring

γ < µ H µ C µ C a HH d H -µ H a CH d C µ H a CC d C -µ C a HC d H (44)
(assuming first µ C a HH d H > µ H a CH d C , a CC µ H d C > a HC µ C d H ), µ H , µ C > 0, (45)
α H µ C r H < α C µ H r C , α H µ H r C < α C µ C r H (46)
(α C r H µ C -α H r C µ H ) (u max 2 ) 2 + 2(r H µ C -r C µ H )u max 2 + α H r H µ C -α C r C µ H α H α C < 0. ( 47
)
Note that the two last assumptions are satisfied as soon as α H α C is very small, at least compared to µ H µ C . This amounts to saying that cytostatic drugs specifically target the cancer cells better than cytotoxic drugs do. This last necessary condition motivates the definitions

φ 1 := µ H p H ρ H + µ C p C ρ C ,
and (abusively, since this quantity also depends on t)

ψ(u 2 ) := r H p H ρ H 1 + α H u 2 + r C p C ρ C 1 + α C u 2 .
Let us first analyse a constrained arc on [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF], whenever it is not reduced to a singleton.

Arc on the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF].

First note that ρ C = ρ C ρ H ρ H = ρ C ρ H θ H ρ 0
H is bounded from above by γ θ H ρ 0 H . If we differentiate the constraint, we find that u 1 and u 2 are determined by

r H 1 + α H u 2 -d H (a HH θ H ρ 0 H + a HC ρ C ) -u 1 µ H = 0, (48)
together with the fact that

(u 1 , u 2 ) ∈ arg max r H p H ρ H 1 + α H u 2 + r C p C ρ C 1 + α C u 2 -(µ H p H ρ H + µ C p C ρ C ) u 1 = arg max p C ρ C µ H r C µ H 1 + α C u 2 - r H µ C 1 + α H u 2
One can check that ( 46) and ( 47) are sufficient conditions to have decrease of the function

u 2 → r C µ H 1+α C u 2 -r H µ C 1+α H u 2 . In particular, u 2 = u max 2 if p C < 0, u 2 = 0 if p C > 0.
Thus, the maximisation condition is equivalent to maximising -φ 1 u 1 if p C does not vanish on the arc. Hence, φ 1 = 0 when this condition on p C is fulfilled. We also obtain u 1 in feedback form along the arc, and when p C does not vanish it is given by:

u b,v 1 := 1 µ H r H 1 + α H v -d H (a HH θ H ρ 0 H + a HC ρ C ) where v = 0 or v = u max 2
depending on the sign of p C . We assume that this is an admissible control, i.e., that it satisfies

0 < u b,v 1 < u max 1 (49)
for v = 0 and v = u max 2 , and any 0 ρ C γθ H ρ 0 H . If p C > 0 and u 2 = 0, the dynamics of ρ C on the arc [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] are given by

dρ C dt = 1 µ H (r b -d b ρ C ) ρ C (50)
with

r d := (r C µ H -r H µ C ) + (a HH d H µ C -µ H a CH d C ) θ H ρ 0 H , d b := (a CC µ H d C -a HC µ C d H ) ,
which we assume to be positive. This autonomous ODE leads to a monotonic behaviour of ρ C . In order to ensure that the boundary control u 1 = u b,0 1 is not enough to prevent the increase of ρ C we assume the following

γ θ H ρ 0 H < r d b d . (51)
The previous hypothesis implies that ρ C will increase on an arc on (11) when p C > 0.

Arc on the constraint [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF].

If we differentiate the constraint, we find that R H = R C , i.e., u 1 and u 2 are related to one another by

r H 1 + α H u 2 -d H ρ H (a HH + γa HC ) -u 1 µ H = r C 1 + α C u 2 -d C ρ H (γa CC + a CH ) -u 1 µ C .
We are now set to prove the result:

Proposition 2. Assume ( 8), ( 32), ( 34), ( 42), ( 44), ( 45), ( 46), ( 47), ( 49), [START_REF] Perthame | Transport equations in biology[END_REF] and that there exists an optimal solution which is the concatenation of free and constrained arcs (either on the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] or ( 10)), with associated times

(t i ) 1 i M .
Then, the last three possible arcs are:

• a boundary arc along the constraint [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF].

• a free arc with controls

u 1 = u max 1 and u 2 = u max 2 ,
• a boundary arc along the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] with

u 2 = u max 2 .
The proof is technical and can be found in Appendix C.

Solution of (OCP) in B T for large T : proof of Theorem 2

Recall that we want to solve (OCP) for controls (u 1 , u 2 ) ∈ B T for large T and small T M 2 , a choice motivated by the previous results. For a given T , we denote

ū(T ) 1 , ū(T ) 2
a choice of optimal values for the constant controls during the first phase.

Theorem 2. Assume the hypotheses of Proposition 2. Then asymptotically in T and for T M 2 small, there exists at least one solution to (OCP) in B T . More precisely, there exists

ūopt 1 , ūopt 2 , T ode 2 , u ode 1 , u ode 2 ∈ BV 0, T ode 2 such that if we define the control (u 1 , u 2 ) by (u 1 , u 2 )(t) = ūopt 1 , ūopt 2 on 0, T -T ode 2 , u ode 1 (t -T + T ode 2 ), u ode 2 (t -T + T ode 2 )
on T -T ode 2 , T then up to a subsequence we have lim

T →+∞ C T (u 1 , u 2 ) - inf (u 1 ,u 2 )∈B T C T (u 1 , u 2 ) = 0,
meaning that (u 1 , u 2 ) is quasi-optimal if T is large enough. Furthermore, on T -T ode 2 , T the optimal trajectory trajectory obtained with (u 1 , u 2 ) is the concatenation of at most three arcs:

• a quasi-boundary arc along the constraint (10),

• a free arc with controls u 1 = u max

1 and u 2 = u max 2 ,
• a quasi-boundary arc along the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF], with u 2 = u max 2 . Remark 8. By quasi-boundary arc, we mean that the quasi-optimal control is such that (ρ H , ρ C ) almost saturates the constraints, i.e., up to an error vanishing as T goes to +∞.

Proof.

Up to a subsequence, still denoted T , we can find ūopt 1 , ūopt

T →+∞ sup [T -T 2 ,T ] max (|ρ H -ρH |, |ρ C -ρC |) = 0, (52)
with the notations of Lemma 5: ρ is obtained from the IDE system, while ρ is obtained from the ODE concentrated on x opt H , x opt C . This is a consequence of a slight refinement of Lemma 5. Indeed, for T large, the IDE is almost concentrated on some x

(T ) H , x (T ) C associated to ū(T ) 1 , ū(T ) 2
. The formulae for these quantities given by Theorem 1 show that

x (T ) H , x (T ) C converges to x opt H , x opt C
, hence the concentration of the IDE on x opt H , x opt C and the result [START_REF] Pisco | Non-Darwinian dynamics in therapy-induced cancer drug resistance[END_REF].

As a consequence, the optimal strategy for the ODE, obtained by Proposition 2 is also optimal for the IDE, up to an error vanishing as T goes to infinity. We denote

T ode 2 T M 2 , u ode 1 , u ode 2 ∈ BV 0, T ode 2
the solutions of this optimal control problem. The last statements of the theorem are then a direct consequence of Proposition 2 and the assumption that T M 2 is small, since the IDE and ODE trajectories are arbitrarily close.

Numerical simulations

In this section, we solve (OCP) numerically in the full class A T . We will compare the results with the previous section, and check that alternative strategies to the one given in Theorem 2 are indeed sub-optimal when T is large.

Numerical simulations of the solution to (OCP)

For a survey on numerical methods in optimal control of ODEs, we refer to [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF].

Here, we use direct methods which consist in discretising the whole problem and reducing it to a "standard" constrained optimisation problem. The IDE structure is dealt with a discretisation in phenotype, which adds to the discretisation in time. The dimension of the resulting optimisation problem becomes larger as the discretisation becomes finer. This method is hence computationally demanding and its numerical implementation requires some care. It relies on combining automatic differentiation and the modelling language AMPL (see [START_REF] Fourer | A modeling language for mathematical programming[END_REF]) with the expert optimisation routine IpOpt (see [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF]). Several different numerical tricks (hot start, numerical refinement, etc) were also needed.

For the simulations, we take θ HC = 0. These simulations clearly indicate that for the chosen numerical data, if T is large enough, then the optimal controls are such that:

• the optimal control u 1 is first equal to 0 on a long arc. Then, on a short-time arc, u 1 = u max 1 and then to a value such that the constraint [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF] saturates;

• the optimal control u 2 has a three-part structure, with a long-time starting arc which is a boundary arc, that is, an arc along which the state constraint ( 10) is (very quickly) saturated. It corresponds to an almost constant value for the control u 2 . The last short-time arc coincides with that of u 1 , and along this arc u 2 = u max .

We denote by t s (T ) the switching time, defined by largest time such that u 1 (t) = 0 for all t < t s (T ).

According to the numerical simulations, as T tends to +∞, both x → n C (t s (T ), x), x → n H (t s (T ), x) converge to (weighted) Dirac masses. Since the controls u 1 and u 2 are almost constant on (0, t s (T )), this is in accordance with Theorem 1. The cancer cell population is then concentrated on a phenotype on which the drugs are very efficient.

More precisely, as T tends to +∞, the optimal strategy seems tends to a two-piece trajectory, consisting of:

• a first long-time arc, along the boundary

ρ H (t)
ρ H (t)+ρ C (t) = δ CH , with u 1 (t) = 0 and with a constant control u 2 , at the end of which the populations of healthy and of cancer cells have concentrated on some given sensitive phenotype;

• a second short-time arc along which the populations of healthy and cancer cells are very quickly decreasing.

We also find that the mapping T -→ ρ C (T ) (where ρ C (T ) is the value obtained by solving (OCP) on [0, T ]) is decreasing. This is because our parameters are such that, once concentrated on a sensitive phenotype, the cancer cell population satisfies a controlled ODE for which there exists a strategy letting ρ C converge to 0. Because our model is exponential, we cannot reach 0 exactly but for very small values of ρ C , one can consider that the tumour has been eradicated.

Remark 9. In order to avoid additional lengthy hypotheses, we did not give conditions under which the strategy established in Theorem 2 can further be identified. However, the numerical solutions show that, for generic parameters, it can be expected that:

• the constant controls on the first phase are such that at the end of the first phase, we have saturation of (10),

• the second phase is of time duration T M 2 and starts with a constrained arc along [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF].

Comparison with clinical settings

As explained before, our results advocate for a first long phase which must be all the more long for an initially heterogeneous tumour (with respect to resistance). They also apply to 'born to be bad' tumours [START_REF] Sottoriva | A big bang model of human colorectal tumor growth[END_REF], with high initial heterogeneity with respect to genes or phenotypes in general. Indeed, the heterogeneity or homogeneity we address here is related to one phenotype defined by resistance towards one category of cytotoxic drug. In this sense, our use of the term heterogeneity is unambiguous, functionally defined, and cannot be superimposed on other more classical uses, defined by the accumulation of mutations, such as in [START_REF] Ding | Clonal evolution in relapsed acute myeloid leukaemia revealed by whole genome sequencing[END_REF][START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF][START_REF] Sottoriva | A big bang model of human colorectal tumor growth[END_REF].

This being said, we are ultimately concerned with the application of our optimal control methods to the improvement of classical therapeutic regimens in which repeated courses of chemotherapy are delivered to patients with cancer. To this end, we keep the previous parameters, that are in particular relevant to represent an initially heterogeneous tumour, and we propose for possible implementation in the clinic a quasi-periodic strategy such as in the example defined below:

• As long as ρ H ρ H +ρ C
θ HC , we follow the drug-holiday strategy by choosing u 1 = ū1 = 0, u 2 = ū2 = 0.5 obtained in the previous numerical simulations.

• Then, as long as ρ H > θ H ρ H (0), we use the maximal amount of drugs. As soon as ρ H = θ H ρ H (0), go back to the drug-holiday strategy.

The implementation is straightforward, Figure 6 shows an example for T = 60. This strategy allows to maintain the tumour size below some upper value and to prevent resistant cells from taking over the whole population. However, the tumour is not eradicated and this strategy is far from being optimal: ρ C (T ) is slightly below 1, to be compared to the value obtained with T = 60 (see Figure5) with the optimal strategy, which is around 1.10 -5 . It is another proof of the importance of a long first phase. It also shows that, at least with our parameters, the last arc on the constraint (11) obtained in the previous simulations is instrumental in view of significantly decreasing the tumour size. To assess the importance of the saturation of the constraint ρ H = θ H ρ H (0), we complement the previous strategy with an arc on this constraint, with u 2 = u max 2 , and adequately chosen feedback control u 1 obtained from the equality dρ H dt = 0. We go back to the drug-holiday strategy as soon as ρ C starts increasing again, since it is a sign that the tumour has become too resistant. We choose T = 100 to have enough cycles; the corresponding results are reported on Figure 7 below. They tend to show that ρ C can be brought arbitrarily close to 0 after enough cycles, meaning that there is a chance for total eradication of the tumour.

Conclusion

Summary of the results

By analysing a controlled integro-differential system of cancer and healthy cells structured by a resistance phenotype, we have mathematically investigated the effect of combined chemotherapeutic (cytotoxic and cytostatic) drugs on a tumour. Since we chose a biologically grounded modelling for the resistance phenomenon and took the healthy tissue into account, our approach is tailored for understanding and circumventing the two main pitfalls in cancer therapy: resistance to drugs and toxicity to healthy tissue. The goal of our analysis was indeed twofold: check that our model can reproduce the possible deleterious effect of chemotherapy when MTDs are used (the standard clinical strategy), and propose alternative (optimised) infusion protocols.

Since MTD can first strongly reduce the size of the tumour which then starts growing again, we addressed the first question through an asymptotic analysis of the model. This was performed in Theorem 1, which showed that both populations converge, while the cells concentrate on some phenotypes. This theorem extends results that so far were established only in the case of a single integro-differential equation, although the ideas are radically different because the usual technique (integration w.r.t. the phenotype to show convergence, and then relying on the exponential nature of convergence, to show concentration) does not work in our context. The proof of convergence and concentration, presented in Section 2, relies instead on a suitably defined Lyapunov function, whose analysis gives the speed of convergence and concentration. Interestingly, the approach could incorporate controls which are not only constant, but also asymptotically constant.

The rest of the article was then devoted to addressing the second question, by considering the optimal control problem (OCP) of minimising the number of cancer cells on a given time interval [0, T ], keeping the tumour size in check and limiting damage to the healthy tissue. In Section 3, we gave several rigorous mathematical arguments to explain why, when T is large, a good strategy is to first steer the cancer cell population on an appropriate phenotype by first giving constant doses for a long time. These arguments justified a restriction to a smaller class of controls for which we managed to identify an asymptotically optimal strategy in large time, presented in Theorem 2.

In Section 4, we showed through numerical simulations that, when T increases, the optimal solution is indeed increasingly close to a two-phase trajectory. The first very long phase consists in giving low doses of drugs in order to let the cancer cell population concentrate on a given sensitive phenotype. The doses are chosen as low as the constraint on the relative tumour size allows it. Our results advocate for a first long phase which must be all the more long for an initially heterogeneous tumour (with respect to resistance). During the second phase, we numerically recover the expected trajectory, given by Theorem 2: high doses are given (MTD as long as the constraint on the healthy tissue does not saturate) and the cancer cell population quickly decreases.

Possible generalisations

We have focused on a 1 dimensional phenotype x ∈ [0, 1]. In applications, however, it might be suitable to consider multi-dimensional phenotypes in order to account for the complexity of resistance. This is for example what is done in [START_REF] Chisholm | Emergence of Drug Tolerance in Cancer Cell Populations: An Evolutionary Outcome of Selection, Nongenetic Instability, and Stress-Induced Adaptation[END_REF] where the relevant phenotype is 2 dimensional and decided to be the combination of proliferation and survival potentials. With some technical adaptations, the results established in this paper generalise to any compact subset of R d , d 1. For the application one would need to specify how the functions depend (monotonically or not) on the various components of the phenotype.

A possible generalisation for our model is to take into account the fact that cells can change phenotype, for instance through (random) genetic mutations, i.e., modifications of the DNA. These are irreversible and are passed from one cell to its daughter cells through division. However, it is now widely believed that such mutations are very rare with respect to the typical timescales that are of medical interest (which in the case of drug resistance phenomena are much shorter than the timescale of a human life, a time during which mutations certainly exist and can explain the development of diseases, see, e.g. [START_REF] Hirsch | Genetic hierarchy and temporal variegation in the clonal history of acute myeloid leukaemia[END_REF] about acute myeloid leukaemia), and thus they can be disregarded at least as a first approximation. In contrast, epimutations (i.e., mechanisms which alter gene expression but not the DNA sequence base pairs themselves) are thought to be much more frequent [START_REF] Brown | Poised epigenetic states and acquired drug resistance in cancer[END_REF][START_REF] Shen | Interplay between the cancer genome and epigenome[END_REF][START_REF] You | Cancer genetics and epigenetics: Two sides of the same coin?[END_REF].

If epimutations are exclusively random, they can be modelled by a diffusion term, in which case (4) would be modified as follows:

∂n H ∂t (t, x) = R H (x, ρ H (t), ρ C (t), u 1 (t), u 2 (t)) n H (t, x) + β H ∂ 2 n H ∂x 2 (t, x), ∂n C ∂t (t, x) = R C (x, ρ C (t), ρ H (t), u 1 (t), u 2 (t)) n C (t, x) + β C ∂ 2 n C ∂x 2 (t, x), (53)
together with Neumann boudary conditions in x = 0 and x = 1. Here, β H and β C stand for the random epimutations rates of the healthy and cancer cell populations, respectively.

The Darwinian idea that the fittest individuals are selected exclusively because of random events affecting the genome or its expression has been recently challenged as observations on genomic evolution cannot be accounted for by sheer Darwinian mechanisms [START_REF] Ling | Extremely high genetic diversity in a single tumor points to prevalence of non-darwinian cell evolution[END_REF], maybe also as ideas from Lamarck are regaining popularity. Such theories advocate the existence of adaptive behaviours: individuals actively adapt to their environment, seeking for phenotype changes that make them fitter. These can be seen as stress-induced epimutations and can be mathematically modelled by an advection term, as in [START_REF] Chisholm | Emergence of Drug Tolerance in Cancer Cell Populations: An Evolutionary Outcome of Selection, Nongenetic Instability, and Stress-Induced Adaptation[END_REF][START_REF] Lorenzi | Effects of an advection term in nonlocal Lotka-Volterra equations[END_REF]

together with Neumann boudary conditions in x = 0 and x = 1. Here, v H and v C are the velocities with which healthy and cancer cells adapt to their environment, respectively, which are assumed to vanish in x = 0 and x = 1. Because we assume that the adaptation of cells is induced by the stress created by the drugs, v H and v C depend on u 1 and u 2 .

Another extension could involve a mixed deterministic/stochastic framework, namely using a piecewise deterministic Markov process (PDMP [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF], see [START_REF] Renault | Optimal control of infinite-dimensional piecewise deterministic markov processes and application to the control of neuronal dynamics via optogenetics[END_REF] for the optimal control of this class of equations). In these models, mutations are stochastic jumps between deterministic (and phenotypically reversible) models, each jump becoming less and less rare in the course of phenotypic evolution in the deterministic processes. Furthermore, in these models, the probability of jump would depend exclusively on (and as an increasing function of) the phenotype structure variable, that would thus bear a quantitative meaning of malignancy, or phenotype plasticity entraining genetic instability (this last point is discussed with references in [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF]).

A final extension should stem from the fact that tumours are also very heterogeneous in space (for example, because cells at the outer rim and cells at the centre of a tumour spheroid encompass very different metabolic conditions; more genenerally, high heterogeneity depending on space has been experimentally shown in solid tumours [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF][START_REF] Sottoriva | A big bang model of human colorectal tumor growth[END_REF], which should lead to also structure the populations of cells according to an added space variable. Another modelling advantage of such representation is that the interaction between the tumour and the healthy tissue is itself spatial, since part of it essentially happens at the boundary of the tumour, through direct contact. For possible cancer models taking both phenotype and space into account, we refer to [START_REF] Jabin | Selection-Mutation dynamics with spatial dependence[END_REF][START_REF] Lorz | Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors[END_REF][START_REF] Mirrahimi | Asymptotic analysis of a selection model with space[END_REF].

Open problems

Asymptotic analysis. In this paper, we have extended well-known results for a single IDE to systems of IDEs. However, in applications it could be interesting to consider the case of general competitive of the form

∂n H ∂t (t, x) = R H (x, ρ H , ρ C ) n H (t, x), ∂n C ∂t (t, x) = R C (x, ρ C , ρ H ) n C (t, x), (55)
with R H decreasing in ρ C , R C decreasing with ρ H . Proving convergence and concentration for such systems is completely open. Indeed, our Lyapunov function is specifically suited to the specific linear setting of the model ( 4): it cannot be applied to any general competitive system. Note that some numerical simulations indicate that no oscillations occur, which may mean that ρ H and ρ C converge.

Similarly, characterising the asymptotic behaviour of PDE systems like ( 53) and (54) (even without the control terms) is an interesting and open problem, even for a single equation. Let us mention that when the rate of mutations is small, much can be found in the literature on the asymptotics of these models when this small parameter goes to 0, after a proper rescaling of time [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF][START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF]. If this parameter is fixed, classical asymptotics for t going to infinity have up to our knowledge not been carried out.

Optimal control. In order to try and solve the optimal control problem (OCP), we had to treat the question of the optimal control of IDEs. Although a PMP exists for such equations, we took another path because the resulting equations were too intricate.

A key idea to justify the restriction to a class of controls which are first constant on a long phase is to prove that the tumour (among all possible tumours of given size) which can be treated the most efficiently is homogeneous in phenotype, i.e., a Dirac mass in mathematical terms: the first long phase then aims at approaching this 'ideal' situation to start the second phase. We established the optimality of Dirac masses for a short time, but not for any time.

If we want to analyse similar optimal control problems for [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'What does not kill me strengthens me[END_REF] or [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF], we shall have to deal with optimal control of PDEs for which techniques are very different [START_REF] Coron | Control and nonlinearity[END_REF][START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]. A first approach would be to focus on [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'What does not kill me strengthens me[END_REF], and see whether and how the optimal controls converge to the ones obtained in this paper as the rates of epimutations go to 0.

A Proof of Lemma 1

Proof. We are going to prove that ρ C is a BV function. To that end, let us prove that ρ C is bounded from above, and that it has integrable negative part.

First step: upper bound for ρ C .

The existence of such a bound comes from integrating the equation with respect to x:

ρ C (t) = 1 0 R C (x, ρ C , 0ū 1 , ū2 ) n C (t, x) dx.
If ρ C is too large, the right hand side is negative, forcing ρ C to decrease. It proves the claim on the upper bound for ρ C .

Similarly, because of assumption ( 14), ρ C increases if ρ C is too close to 0: ρ C is bounded from below by some ρ min C > 0.

Second step: estimate on the negative part of ρ C .

We define q C := dρ C dt and wish to prove that (q C ) -∈ L 1 (0, +∞). We differentiate

dρ C dt = 1 0 n C R C to obtain: dq C dt = 1 0 n C R 2 C + 1 0 n C ∂R C ∂ρ C q C
It provides an upper bound for the negative part of q C :

d(q C ) - dt 1 0 n C ∂R C ∂ρ C (q C ) --a CC d min C ρ min C (q C ) -
where 0 < d min C d C on [0, 1]. We conclude that the negative part of q C vanishes exponentially (and consequently, is integrable over the half-line). Therefore, ρ C converges to some ρ ∞ C > 0.

Figure 1 :

 1 Figure 1: Simulation with u 1 (t) = Cst = 3.5 and u 2 (t) = Cst = 2, in time T = 10. At the top, left and middle: evolution in time of the curves x → n H (t, x) and x → n C (t, x), with the initial conditions in black, and the final ones in red. At the right, top and bottom: graphs of t → ρ C (t) and of t → ρ H (t). At the bottom, left and middle: graphs of t → ρ H (t)

Figure 2 :

 2 Figure 2: Former function µ C in blue, and new function µ C in red.

Figure 3 : 1 , u max 2 )

 312 Figure 3: Simulation with u 1 (t) = Cst = 3.5 and u 2 (t) = Cst = 2, in time T = 10, with the new function µ C .

2 as

 2 T → +∞. These values for the constant controls yield asymptotic phenotypes x opt H , x opt C thanks to Theorem 1. Then, for any choice of time T 2 T M 2 and BV controls (u 1 , u 2 ) on (T -T 2 , T ), lim

  [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF], θ H = 0.6, = 0.1. We let T take the values T = 30 and T = 60. The results are reported on Figures4 and 5respectively.

Figure 4 :

 4 Figure 4: Simulation of (OCP) for T = 30.

Figure 5 :

 5 Figure 5: Simulation of (OCP) for T = 60.

2

 2

Figure 6 :

 6 Figure 6: Quasi-periodic strategy, for T = 60.

Figure 7 :

 7 Figure 7: Second quasi-periodic strategy, for T = 100.

  Third step: identification of ρ ∞ C . Now, we have n C (t, x) = n 0 C (x) exp r C (x) 1 + α C ū2 -ū1 µ C (x) t -d C (x)a CC t 0 ρ C (s) dsFor large t, we have t 0 ρ C (s) ds ∼ ρ ∞ C t and hence the asymptotic behaviour depends on the function b C defined on [0, 1] by bC (x) = r C (x) 1 + α C ū2 -ū1 µ C (x) -d C (x)a CC ρ ∞ C .

  . This would lead to a model of the form ∂n

	∂n C ∂t	(t, x) +	∂ ∂x	(v

H ∂t (t, x) + ∂ ∂x (v H (x, u 1 (t), u 2 (t)) n H (t, x)) = R H (x, ρ H (t), ρ C (t), u 1 (t), u 2 (t)) n H (t, x) + β H ∂ 2 n H ∂x 2 (t, x), C (x, u 1 (t), u 2 (t)) n C (t, x)) = R C (x, ρ C (t), ρ H (t), u 1 (t), u 2 (t)) n C (t, x) + β C ∂ 2 n C ∂x 2 (t, x),

A boundary arc (for the state constraint g (ρH , ρC ) 0) is an arc along which g (ρH , ρC ) = 0, i.e., the constraint is saturated. A free arc is an arc along which g (ρH , ρC ) < 0, i.e., the constraint is not saturated.

The starting and ending points of a boundary arc are called junction points if they are distinct, and contact points if they coincide (i.e., if the arc is reduced to a singleton).

Let B C ⊂ [0, 1] be the set of points at which the function b C reaches its maximum.

Let us prove that b C (x) = 0, for every x ∈ B C . We argue by contradiction. If b C (x) > 0 for some x ∈ B C , then there exists a nontrivial interval I ⊂ [0, 1] along which b C is positive, and therefore n C (t, •) takes larger and larger values along I as t increases. This contradicts the fact that ρ C (t) converges to ρ ∞ C . Similarly, if b C < 0 globally, ρ C converges to 0, a contradiction.

The function b C is thus nonpositive on [0, 1], and vanishes at any point of B C . The lemma follows easily.

B Proofs for the simplified optimal control problems B.1 Proof of Lemma 3 Proof. Using the family u defined in Remark 5, we obtain the corresponding ρ (T ), which can be computed exactly, as well as its limit. It is given by ρ opt (T ) := ρ opt (T -) exp(-µu 1,max ) where ρ opt is the function obtained through d dt ρ opt (t) = (r -dρ opt (t))ρ opt (t) for t < T , and ρ opt (0) = ρ 0 . Now, let any u satisfy [START_REF] Greaves | Evolutionary determinants of cancer[END_REF]. The solution of [START_REF] Greaves | Cancer stem cells: back to Darwin?[END_REF] with u is thus a subsolution of that satisfied by ρ opt , leading to ρ ρ opt on [0, T ). Using u 0, we also have

Since ρ opt (T -) = ρ 0 exp T 0 (r -dρ opt (s)) ds and ρ ρ opt , this implies ρ(T ) ρ opt (T ). Let us now investigate the possible case of equality to prove that the infimum is not attained: the foregoing equality implies that there is equality if and only if T 0 u ds = u 1,max (the contraint is saturated) and exp T 0 (r -dρ(s)) ds = exp T 0 (r -dρ opt (s)) ds , whence ρ ≡ ρ opt on [0, T ). As ρ is continuous, ρ(T ) would be given by taking u ≡ 0, which is not optimal.

B.2 Proof of Lemma 4

Proof. To account for the L 1 constraint (C1), we augment the system by defining another state variable y, whose dynamics are given by dy dt = u, leading to:

The constraint (C1) thus rewrites y(T ) u 1,max .

According to the Pontryagin maximum principle (see [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]), there exist absolutely continuous adjoint variables p ρ and p y on [0, T ], and p 0 0, such that:

where the Hamiltonian is

Thus, p y is some constant, and p ρ does not change sign on [0, T ].

The maximisation of the Hamiltonian leads to defining the switching function φ := p yµp ρ ρ. u is thus equal to u ∞,max whenever φ > 0, equal to 0 whenever φ < 0.

The transversality condition is that the vector p ρ p y (T ) -p 0 1 0 must be orthogonal to the tangent space of {(p, y) ∈ R 2 , y u 1,max } at the point (ρ(T ), y(T )).

First case.

If y(T ) < u 1,max , then the transversality conditions imply p ρ (T ) = p 0 and p y ≡ 0. p 0 = 0 since otherwise we would have (p ρ , p y , p 0 ) ≡ 0. Thus, in this case, p ρ (T ) < 0 and p ρ is negative on the interval [0, T ]. The switching function φ is therefore positive on the whole [0, T ], which would imply u ≡ u ∞,max . This is a contradiction since a consequence is

If y(T ) = u 1,max , we still have p ρ (T ) = p 0 . As in the first case, we cannot have p y = 0.

Let us first remark that φ cannot be positive nor negative on the whole interval, since otherwise u ≡ u ∞,max , a contradiction, or u ≡ 0, which is clearly not optimal. If p 0 = 0, p ρ ≡ 0, so that φ has the sign of p y = 0, a contradiction. Therefore, p ρ < 0 on [0, T ] as before, and this implies p y < 0 to ensure that φ changes sign.

The derivative of φ is given by dφ dt = -µdp ρ ρ 2 > 0. Thus, φ is increasing and u is bangbang with one switching only. The fact that y(T ) = T 0 u(s) ds = u 1,max imposes that this switching happens at T 1 (T ) as announced, which ends the proof.

C Proof of Proposition 2

Proof. If the constraint [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF] does not saturate on the whole [0, t f ], we distinguish on whether the last arc is a free arc or a boundary arc on [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF].

First case: the last arc is a boundary arc on [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF], not reduced to a singleton.

In this case, t M = t f and there can be a jump on the adjoint vector at t f . Let us start by proving the following: Lemma 6. p 0 < 0.

Proof. We argue by contradiction and assume p 0 = 0. We first look at the interval [t M -1 , t f ], and assume, also by contradiction, that ν M > 0. Then p

From assumption [START_REF] Perthame | Transport equations in biology[END_REF], this means that ρ C decreases locally around t f , a contradiction since t f is free (a better strategy would be to stop before ρ C starts increasing): ν M = 0. Now, let us prove that p H , p C and η 1 vanish identically on

we define the maximal interval [t 0 , t ) on which p C > 0 (resp., p C < 0), with p C (t ) = 0. In this case, we know that the switching function φ 1 vanishes on [t 0 , t ], hence p H factorises with p C . Coming back to the equation on p C , we have p C = β C p C on (t 0 , t ), for some function β C . Since p C (t ) = 0, this imposes p C ≡ 0 on the interval, a contradiction. Thus p C is identically 0 on the whole (t M -1 , t f ), and so are p H (from the equation on p C ) and η 1 (from the equation on p H ).

We now analyse the arc [t M -2 , t M -1 ]. From the previous step, we know that φ 1 (t M -1 ) = 0.

. By continuity of the Hamiltonian, we get R H t - M -1 = 0. At the left of t M -1 , u 1 and u 2 saturate at their maximal values. At the right of t M -1 , R H = 0 but this imposes u 1 < u max 1 or u 2 < u max 2 since, owing to [START_REF] Jabin | On selection dynamics for competitive interactions[END_REF], ρ H decreases for the maximal values. Thus, 0 = R H t - M -1 < R H (t M -1 ) = 0, a contradiction. Finally, we have proved ν M -1 = 0.

Standard Cauchy-Lispchitz arguments, together with the result p H (t - M -1 ) = p C (t - M -1 ) = 0 yield that p H and p C are also identically null on the interval [t M -2 , t M -1 ]. Repeating these arguments on the whole [0, t f ], we find that p, p 0 , η 1 , η 2 and the (ν i ) i=1,...,M are all zero, in contradiction with condition 1 given by the PMP (see Section 3). Thus p 0 < 0 and we set p 0 = -1. This normalisation is allowed because the final adjoint vector (p(t f ), p 0 ) is defined up to scaling. Again, we start by analysing the PMP on

and φ 1 = 0 locally around t f . This implies p H > 0 also locally around t f . In particular, ν M = 0. Using the same reasoning as before with p C = β C p C , we get this time that p C and p H have constant sign on (t M -1 , t f ): p C < 0 and p H > 0.

Let us now first assume ν

. At t M -1 , we would thus have continuity of u 2 and not u 1 since u 1 < u max 1 on [t M -1 , t f ] from assumption [START_REF] Navin | Inferring tumor progression from genomic heterogeneity[END_REF]. In such a case, it holds true that there can be no jump on the adjoint vector (see for instance [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF]), which is contradictory unless ν M -1 is such that p H t - M -1 > 0, which we assume from now on.

Let us now analyse the interval [t M -2 , t M -1 ], on which we will prove that

Because η 1 and η 2 vanish on such an interval, it is easy to prove from standard Cauchy-Lipschitz uniqueness arguments that p C < 0 and p H > 0 on [t M -2 , t M -1 ]. Also, because of (44) the inequality

is satisfied on [0, t f ]. Let us prove that this implies φ 1 < 0 on (t M -2 , t M -1 ). For that purpose, we will prove that whenever φ 1 (t 0 ) = 0, its derivative satisfies φ 1 (t 0 ) > 0. Note that we already know that φ 1 (t - M -1 ) φ 1 (t M -1 ) = 0. For such a time t 0 we indeed obtain

Combined with [START_REF] Scharovsky | Metronomic chemotherapy: changing the paradigm that more is better[END_REF], this yields φ 1 (t 0 ) > 0, as announced. Thus

For u 2 , the proof is a bit more involved because the dependence is not linear. In what follows, we generically denote

as a function of u 2 , whose derivative has the opposite sign of P (u 2 ), where

0 on (0, t f ). We consider two cases, depending on the sign of φ (α C r H , α H r C ) . Note that (8) implies the order [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF] and φ 1 < 0, P (0

Let us first assume φ (α C r H , α H r C ) < 0, in which case all the coefficients of the polynomial are negative. Let us denote u + the greater root of this polynomial. Since the coefficient in front of u 2 is negative, the function ψ is increasing with u 2 on (u + , +∞). We cannot have u + 0 because of the signs of the coefficients: u max 2 maximises the function of interest. If φ (α C r H , α H r C ) = 0, it is easy to see that the same result holds.

Now, let us assume that φ

), we find

We conclude that P (u max

2

) < 0 thanks to [START_REF] Lorz | Long-term behaviour of phenotypically structured models[END_REF].

Thus, we have proved that, on (t M -2 , t M -1 ), u 1 = u max 1 and u 2 = u max

2

. Note that the result actually implies ν M -1 = 0. However the same reasoning with ν M -1 = 0 works and we obtain u 1 = u max 1 and u 2 = u max 2 . From assumption (34), ρ C ρ H increases backwards. If this ratio reaches the value γ, i.e., if the system saturates the constraint (10) (if not, t M -2 = 0), then we have a potential boundary arc on [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF] on (t M -3 , t M -2 ).

Second case: the last arc is a boundary arc on [START_REF] Busse | Mass concentration in a nonlocal model of clonal selection[END_REF], reduced to a singleton.

Note that, again, t M = t f . This case is handled as the previous one: p 0 cannot be 0 and φ 1 (t - f ) 0. Because of this result, the whole reasoning made above in the previous case applies: there is an unconstrained arc with u 1 = u max 1 and u max 2 . If there is a previous arc, it is a constrained arc on [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF].

Third case: the last arc is a free arc.

Again, the same kind of arguments are enough to prove that p 0 < 0, and u 1 = u max 1 and u max 2 on this arc. If there is a previous arc, it is a constrained arc on [START_REF] Burrell | The causes and consequences of genetic heterogeneity in cancer evolution[END_REF].