Luís Almeida

Patrizia Bagnerini

Giulia Fabrini

Barry D Hughes

Tommaso Lorenzi

Evolution of cancer cell populations under cytotoxic therapy and treatment optimisation: insight from a phenotype-structured model *

We consider a phenotype-structured model of evolutionary dynamics in a population of cancer cells exposed to the action of a cytotoxic drug. The model consists of a nonlocal parabolic equation governing the evolution of the cell population density function. We develop a novel method for constructing exact solutions to the model equation, which allows for a systematic investigation of the way in which the size and the phenotypic composition of the cell population change in response to variations of the drug dose and other evolutionary parameters. Moreover, we address numerical optimal control for a calibrated version of the model based on biological data from the existing literature, in order to identify the drug delivery schedule that makes it possible to minimise either the population size at the end of the treatment or the average population size during the course of treatment. The results obtained challenge the notion that traditional high-dose therapy represents a 'one-fits-all solution' in anticancer therapy by showing that the continuous administration of a relatively low dose of the cytotoxic drug performs more closely to the optimal dosing regimen to minimise the average size of the cancer cell population during the course of treatment.

Introduction

The emergence of resistance to cytotoxic drugs (i.e. therapeutic agents commonly used to kill cancer cells) is a major obstacle to successful cancer chemotherapy and management of disease relapse [START_REF] Brown | Poised epigenetic states and acquired drug resistance in cancer[END_REF][START_REF] Dagogo-Jack | Tumour heterogeneity and resistance to cancer therapies[END_REF][START_REF] Gottesman | Mechanisms of cancer drug resistance[END_REF][START_REF] Holohan | Cancer drug resistance: an evolving paradigm[END_REF][START_REF] Turner | Genetic heterogeneity and cancer drug resistance[END_REF]. Therefore, understanding the mechanisms which orchestrate the emergence of resistance to cytotoxic therapy is a timely and key challenge for the cancer research community [START_REF] Chabner | Chemotherapy and the war on cancer[END_REF].

A synergistic approach involving empirical and theoretical work has provided evidence in favour of the idea that the emergence of resistance to cytotoxic drugs is, in essence, an evolutionary process driven by the selective pressure that such drugs exert on tumour cells and other cells of the tumour microenvironment [START_REF] Aktipis | Life history trade-offs in cancer evolution[END_REF][START_REF] Dexter | Heterogeneity of tumor cells from a single mouse mammary tumor[END_REF][START_REF] Greaves | Clonal evolution in cancer[END_REF][START_REF] Junttila | Influence of tumour micro-environment heterogeneity on therapeutic response[END_REF][START_REF] Korolev | Turning ecology and evolution against cancer[END_REF][START_REF] Nowell | The clonal evolution of tumor cell populations[END_REF]. For this reason, increasing attention has been given to integrodifferential equations (IDEs) and nonlocal partial differential equations (PDEs) modelling adaptive dynamics in populations structured by physiological traits as possible tools to study in silico the dynamics of cancer cells and their response to cytotoxic drugs [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: biological and mathematical assessment, theoretical treatment optimisation[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]. These models have been proven useful to test biological hypotheses and complement experimental research by enabling extrapolation beyond scenarios which can be investigated through in vitro and in vivo experiments [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF][START_REF] Delitala | A mathematical model for progression and heterogeneity in colorectal cancer dynamics[END_REF][START_REF]A mathematical model for the dynamics of cancer hepatocytes under therapeutic actions[END_REF][START_REF] Hodgkinson | Computational approaches and analysis for a spatio-structural-temporal invasive carcinoma model[END_REF][START_REF] Lavi | The role of cell density and intratumoral heterogeneity in multidrug resistance[END_REF][START_REF]Simplifying the complexity of resistance heterogeneity in metastasis[END_REF][START_REF] Lorenzi | The role of spatial variations of abiotic factors in mediating intratumour phenotypic heterogeneity[END_REF][START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF]. Moreover, the results presented by Olivier & Pouchol [START_REF] Olivier | Combination of direct methods and homotopy in numerical optimal control: Application to the optimization of chemotherapy in cancer[END_REF] and Pouchol et al. [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF] have highlighted how optimal control of such IDE and PDE models can both inform the design of optimised anticancer treatments and raise new interesting mathematical questions. In fact, while the literature on the optimal control of ordinary differential equations (ODEs) modelling the dynamics of cancer cell populations under therapy is now vast -we refer the interested reader to the recent monograph by Schättler & Ledzewicz [START_REF] Schättler | Optimal control for mathematical models of cancer therapies[END_REF] and references therein -there has been little prior work on the optimal control of physiologically structured population models formulated in terms of IDEs and nonlocal PDEs.

Lorenzi et al. [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF] have recently presented a phenotype-structured model of evolutionary dynamics in a population of cancer cells exposed to the action of a cytotoxic drug. In this model, the cell phenotypic state is represented by a continuous variable which is related to the expression level of a gene that controls the cellular level of cytotoxic-drug resistance. The model consists of a nonlocal parabolic equation that governs the evolution of the cell population density function. Phenotypic variations are modelled through a diffusion operator, while a nonlocal reaction term is used to take into account the effect of cell proliferation and death under the selective pressure of the cytotoxic drug. The long-time asymptotic behaviour of the solutions to the model equation have been studied [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF] in the absence of the cytotoxic drug and in the presence of a constant drug dose. The results obtained have provided some insight into the relative contributions of phenotypic variations, intracellular competition for resources and the cytotoxic drug as drivers of adaptation in cancer cell populations. However, three key questions remain open in the mathematical framework of this model.

(i) How do the size and the phenotypic composition of the cell population evolve when the drug dose changes with time?

(ii) What is the optimal drug delivery schedule that allows one to minimise the number of viable cancer cells at the end of the treatment?

(iii) How does the optimal drug delivery schedule change in the case where one wishes to minimise the average size of the cancer cell population during the course of treatment?

In this paper, we address question (i) by developing a novel method for constructing exact solutions to the nonlocal parabolic equation that governs the evolution of the cell population density function. Using this method, we derive a semi-explicit expression for the cell population density function when the administered drug dose is a generic function of time, and we make this expression fully explicit in the case where the drug dose is piecewise-constant in time. This allows us to perform a systematic investigation of the way in which the size and the phenotypic composition of the cancer cell population change in response to variations of the drug dose and other evolutionary parameters.

In order to address questions (ii) and (iii), we consider a calibrated version of the model, which is based on biological data from the existing literature, and we solve numerically different optimal control problems whereby the size of the cancer cell population (i.e. the integral of the population density function) represents the state variable, the delivered drug dose is chosen as control variable and the cost functional is defined by either the population size at the end of the treatment or the average population size during the course of treatment. Such optimal control problems are subject to the dynamic constraints given by the nonlocal parabolic equation of the model and some control constraints which translate into mathematical terms the clinical limitations imposed by the side effects of cytotoxic drugs. The results obtained challenge the notion that traditional high-dose therapy represents a 'one-fits-all solution' in anticancer therapy by showing that the continuous administration of a relatively low dose of the cytotoxic drug performs more closely to the optimal dosing regimen in the case where the goal is to minimise the average size of the cancer cell population during the course of treatment.

The model

We consider a well-mixed population of cancer cells structured by the expression level of a gene which is linked both to the cell proliferation rate and to the cellular level of cytotoxic-drug resistance -such as the genes ALDH1, CD44, CD117 or MDR1 [START_REF] Hanahan | Hallmarks of cancer: the next generation[END_REF][START_REF] Medema | Cancer stem cells: the challenges ahead[END_REF]. Cells inside the population proliferate or die, compete for limited resources, and undergo phenotype variations due to epimutations -i.e. heritable changes in gene expression that leave the sequence of bases in the DNA unaltered [START_REF] Oey | On the meaning of the word 'epimutation[END_REF]. In particular, we focus on the case of spontaneous epimutations, that is, epimutations that occur randomly due to nongenetic instability and are not induced by any selective pressure [START_REF] Huang | Genetic and non-genetic instability in tumor progression: link between the fitness landscape and the epigenetic landscape of cancer cells[END_REF]. Furthermore, a cytotoxic drug can be present, which acts by increasing the death rate of cancer cells.

We represent the expression level of the gene under consideration by a continuous variable y ∈ R ≥0 . Moreover, following the biological ideas proposed by Pisco & Huang [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'what does not kill me strengthens me[END_REF], we assume that there is a level of expression y H which endows cells with the highest level of cytotoxic-drug resistance and a level of expression y L < y H which gives to the cells the highest proliferation rate when the drug is not present. We model the phenotypic state of each cell in the population by means of the rescaled variable x ∈ R with

x = y -y L y H -y L ,
so that the state x = 1 corresponds to the highest level of cytotoxic-drug resistance, whereas the state x = 0 corresponds to the highest proliferation rate in the absence of the drug. At any time t ∈ [0, T ], with T > 0 being a generic final time, we describe the number density of cells in the phenotypic state x at time t (i.e. the cell population density) by means of the function n(x, t) ≥ 0, and we use the function u(t) ≥ 0 to model the (rescaled) dose of the cytotoxic drug. Moreover, we compute the total number of cells (i.e. the population size) ρ(t) and the cell mean phenotypic state µ(t), respectively, as

ρ(t) = R n(x, t) dx and µ(t) = 1 ρ(t) R
x n(x, t) dx.

(2.1)

The evolution of the population density function n(x, t) is governed by the following nonlocal parabolic equation:

∂n ∂t = R x, ρ(t), u(t) n proliferation/death + β ∂ 2 n ∂x 2 , spontaneous epimutations (x, t) ∈ R × [0, T ], (2.2)
which can be obtained from stochastic individual-based models that track the dynamics of single cells via meanfield approximation [START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to population-level consequences[END_REF], through probabilistic methods in the limit of large cell numbers [START_REF] Champagnat | The canonical equation of adaptive dynamics: a mathematical view[END_REF][START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF] or by using the principle of mass conservation [START_REF] Stace | A phenotype-structured individual-based model for the evolution of cancer cell populations under chemotherapy[END_REF].

In equation (2.2), the diffusion term models the effects of spontaneous epimutations, which occur at rate β ∈ R >0 [START_REF] Becker | Spontaneous epigenetic variation in the arabidopsis thaliana methylome[END_REF][START_REF] Navin | Cancer genomics: one cell at a time[END_REF]. The reaction term takes into account the effects of cell proliferation, natural death, competition for resources and the cytotoxic action of the drug. The functional R x, ρ(t), u(t) represents the fitness of cancer cells in the phenotypic state x under the environmental conditions determined by the population size ρ(t) and by the drug concentration u(t). For avoidance of doubt, we emphasis that u(t) relates to the delivered drug dose. In fact, building upon the modelling approach employed in previous papers [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF][START_REF] Lavi | The role of cell density and intratumoral heterogeneity in multidrug resistance[END_REF][START_REF]Simplifying the complexity of resistance heterogeneity in metastasis[END_REF][START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], we do not describe the fine details of clinical administration protocols and how they relate to the kinetics of drug uptake, or the elimination of the drug over time by excretion.

Building upon the modelling strategies presented in [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF], we make use of the following definition for the fitness functional: R x, ρ(t), u(t

) = p(x) -d ρ(t) -k(x, u(t)). (2.3)
In the definition given by equation (2.3), the function p(x) stands for the net proliferation rate of cancer cells (i.e. the difference between the rate of cell division and the rate of natural death) in the phenotypic state x, while the function k(x, u(t)) is the rate of death caused by the cytotoxic drug. Moreover, the saturating term d ρ(t) translates into mathematical terms the idea that higher total numbers of cells correspond to less available resources and space in the system, and thus to more intense intrapopulation competition. The parameter d ∈ R >0 models the rate of cell death due to intrapopulation competition. Based on the ideas proposed by Chisholm et al. [START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to population-level consequences[END_REF] and Lorenzi et al. [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF][START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF], we define the functions p and k as

p(x) = γ -η x 2 , k(x, u) = (x -1) 2 u. (2.4)
In the definitions given by equation (2.4), the parameter γ ∈ R >0 corresponds to the maximum fitness of cancer cells, and the non-linear selection gradient η ∈ R >0 provides a measure of the strength of natural selection in the absence of the cytotoxic drug [i.e. when u(t) = 0]. The definition of the function p(x) is such that if u(t) = 0 the cells in the phenotypic state x = 0 will have the highest fitness. Moreover, the fact that the net proliferation rate p(x) can become negative for values of x sufficiently far from zero captures the idea that phenotypic variants with a low level of fitness cannot survive within the population. The definition of the function k(x, u(t)) models the fact that cells in the phenotypic state x = 1 are fully resistant to the cytotoxic drug and, for cells in phenotypic states other than the most resistant one, the rate of death induced by the drug increases with the drug dose. Besides satisfying these generic properties, the definitions given by equation (2.4) ensure analytical tractability of the model. Moreover, they lead naturally to a smooth fitness functional that is close to the approximate fitness landscapes which can be inferred from experimental data through regression techniques [START_REF] Otwinowski | Inferring fitness landscapes by regression produces biased estimates of epistasis[END_REF].

Remark 1. A little algebra shows that the definitions (2.3) and (2.4) imply

R x, ρ(t), u(t) = γ - η u(t) η + u(t) -(η + u(t)) x - u(t) η + u(t) 2 -d ρ(t).
(2.5)

Therefore, the fittest phenotypic state (i.e. the phenotypic state with the highest fitness) at time t is

x fit (t) = u(t) η + u(t)
and the gene expression level corresponding to the fittest phenotypic state y fit (t) is given by the following equation

y fit (t) = x fit (t) y H + (1 -x fit (t)) y L .
Hence, although y ∈ R ≥0 and x ∈ R, we have that y L ≤ y fit (t) ≤ y H and 0 ≤ x fit (t) ≤ 1 for all t ∈ [0, T ].

Exact solutions

In the case where the administered dose of the cytotoxic drug remains constant over time -i.e. when u(•) = ū with ū ∈ R ≥0 -we have previously shown [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF] that for any initial condition n(x, 0) satisfying the biologically reasonable assumptions

n(x, 0) ∈ L 1 ∩ L ∞ (R), supp (n(•, 0)) ≡ Ω, (3.1)
with Ω ⊂ R being a compact set, the long-time behaviour of the solution n(x, t) to equation (2.2) is given by the following theorem.

Theorem 1. Under assumptions (2.3), (2.4) and (3.1), the integral ρ(t) of the solution to equation (2.2) satisfies the following alternative:

(i) if γ - η ū η + ū -[β(η + ū)] 1/2 ≤ 0 then lim t→∞ ρ(t) = 0; (ii) if γ - η ū η + ū -[β(η + ū)] 1/2 > 0 then lim t→∞ ρ(t) = ρ ∞ ū > 0, with ρ ∞ ū = 1 d γ - η ū η + ū -[β(η + ū)] 1/2 . (3.2) Moreover, in case (ii) equation (2.2) admits a unique nonnegative nontrivial equilibrium solution n ∞ ū (x) with n ∞ ū (x) = ρ ∞ ū √ 2π f ∞ ū β 1/4 exp - 1 2 f ∞ ū β 1/2 (x -µ ∞ ū ) 2 , f ∞ ū = η + ū, µ ∞ ū = ū η + ū . (3.3)
In summary, the asymptotic results established by Theorem 1 formalise the following ideas:

-higher rates of spontaneous epimutations can drive the cancer cell population to extinction under strong selective pressures;

-when there is no cytotoxic drug, the population evolves to be mainly composed of highly proliferative cells;

-higher drug doses correspond to lower cell numbers and promote a selective sweep towards more resistant phenotypic variants.

We refer the interested reader to the work by Lorenzi et al. [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF] for a detailed discussion of the biological implications of these mathematical results.

Remark 2. We shall assume that γ > (β η) This assumption is in accord with biological data from the existing literature (cf. the parameter values listed in Table 1) and captures the fact that, in well-adapted cell populations, the rates of phenotypic variation and the selection gradients are typically smaller than the maximum fitness. If this restriction is not imposed, then there is no choice of ū > 0 for which case (ii) in Theorem 1 applies. Under assumption (3.4), there exists an interval 0 ≤ ū < u † in which case (ii) applies, where u † is the smallest positive real solution of the following algebraic equation for z γ -

η z η + z -[β(η + z)] 1/2 = 0. (3.5)
If we express the dose of cytotoxic drug in terms of the LDα, i.e., the constant value ū of the function u(t) that is required to reduce the equilibrium value of the total number of cells ρ ∞ ū=0 given by equation (3.2) by α%, then u LD100 = u † . In order to take into account toxicity constraints, we shall also assume throughout this paper that

u L ∞ ([0,T ]) < u LD100 . (3.6)
We consider initial conditions of the form

n(x, 0) = ρ 0 √ 2π f 0 β 1/4 exp - 1 2 f 0 β 1/2
x -µ 0 2 with ρ 0 , f 0 ∈ R >0 and µ 0 ∈ R.

This ensures that our discussion includes (but is not limited to) the case in which the system starts in an equilibrium state corresponding to some choice of the parameters β, γ, η and d and either no drug, or some constant drug level below u LD100 . Although the restriction ρ 0 ≤ ρ ∞ ū=0 corresponds to the most relevant biological case, our results remain valid so long as ρ 0 > 0. Furthermore, notice that the factor 1/β can obviously be absorbed into the constant f 0 .

In this section, we construct exact solutions to the Cauchy problem defined by equation (2.2) complemented with the initial condition (3.7). In more detail, we firstly derive a semi-explicit solution in the case where the drug dose is a generic function of time (vid. Proposition 1). Secondly, we make this solution explicit in the case where u(t) is either piecewise-constant (vid. Proposition 2) or constant (vid. Proposition 3).

In the case where the drug dose is a generic function of time, a semi-explicit solution to the nonlocal parabolic equation (2.2) complemented with the initial condition (3.7) is given by the following proposition, the proof of which is provided in Appendix A.

Proposition 1 (Semi-explicit solution for generic u(t)). Under assumptions (2.3), (2.4), (3.4) and (3.7), equation (2.2) admits the exact solution

N (x, t) = ρ(t) √ 2π f (t) β 1/4 exp - 1 2
f (t) β 1/2 (x -µ(t)) 2 , (3.8)
with f (t), µ(t) and ρ(t) satisfying the Cauchy problem

             f (t) = 4(βf (t)) 1/2 (η + u(t) -f (t)), t ∈ (0, T ] µ (t) = 2 β f (t) 1/2 (η + u(t)) u(t) η + u(t) -µ(t) , t ∈ (0, T ] ρ (t) = (Q(t) -dρ(t)) ρ(t), t ∈ (0, T ] ρ(0) = ρ 0 , f (0) = f 0 , µ(0) = µ 0 , (3.9) with Q(t) = γ - η u(t) η + u(t) -(η + u(t)) u(t) η + u(t) -µ(t) 2 -(η + u(t)) β f (t) 1/2 .
(3.10)

Moreover, solving the differential equation for ρ(t) gives

ρ(t) = ρ(0) exp t 0 Q(s) ds 1 + d ρ(0) t 0 exp z 0 Q(s) ds dz -1 . (3.11)
In the case where the drug dose is a piecewise-constant function of time, that is, by the following recursive schemes for i = 0, . . . , N -1 and t ∈ (τ i , τ i+1 ] :

u(0) = u 0 , u(t) = u i 1 (τi,τi+1] (t) for i = 0, . . . , N -1 (3.12) with u i ∈ R ≥0 , τ 0 = 0, τ N = T and τ 0 < τ 1 < . . . < τ N -1 < τ N ,
f (t) =                (η + u i ) tanh 2 2β 1/2 (η + u i ) 1/2 (t -τ i ) + arctanh f (τ i ) η + u i 1/2 , if f (τ i ) < η + u i η + u i , if f (τ i ) = η + u i (η + u i ) coth 2 2β 1/2 (η + u i ) 1/2 (t -τ i ) + arctanh η + u i f (τ i ) 1/2 , if f (τ i ) > η + u i (3.13)
and

µ(t) = u i η + u i + µ(τ i ) - u i η + u i φ i (t), (3.14)
where

φ i (t) =                  ζ i cosech 2β 1/2 (η + u i ) 1/2 (t -τ i ) + arctanh f (τ i ) η + u i 1/2 , if f (τ i ) < η + u i exp -2β 1/2 (η + u i ) 1/2 (t -τ i ) , if f (τ i ) = η + u i λ i sech 2β 1/2 (η + u i ) 1/2 (t -τ i ) + arctanh η + u i f (τ i ) 1/2 , if f (τ i ) > η + u i (3.15)
with

ζ i = f (τ i ) (η + u i ) -f (τ i ) 1/2 and λ i = f (τ i ) f (τ i ) -(η + u i ) 1/2 . ( 3.16)
Remark 3. Under the assumptions of Proposition 2, one can also obtain a recursive scheme for the function Q(t) in equation (3.10) and its integral, which is needed in the computation of ρ(t) via equation (3.11). This is shown by the calculations presented in Appendix D.

In the case where the drug dose is constant over time, that is,

u(t) = ū ∈ R ≥0 for all t ∈ [0, T ], (3.17)
an explicit solution to the nonlocal parabolic equation (2.2) complemented with the initial condition (3.7) is given by the following proposition, the proof of which is given in Appendix C.

f (t) =                (η + ū) tanh 2 2β 1/2 (η + ū) 1/2 t + arctanh f 0 η + ū 1/2 , if f 0 < η + ū η + ū, if f 0 = η + ū (η + ū) coth 2 2β 1/2 (η + ū) 1/2 t + arctanh η + ū f 0 1/2 , if f 0 > η + ū, (3.18) µ(t) =                    ū η + ū + ζ µ 0 - ū η + ū cosech 2β 1/2 (η + ū) 1/2 t + arctanh f 0 η + ū 1/2 , if f 0 < η + ū ū η + ū + µ 0 - ū η + ū exp -2β 1/2 (η + ū) 1/2 t , if f 0 = η + ū ū η + ū + λ µ 0 - ū η + ū sech 2β 1/2 (η + ū) 1/2 t + arctanh η + ū f 0 1/2 , if f 0 > η + ū, (3.19
) where

ζ = f 0 (η + ū) -f 0 1/2 and λ = f 0 f 0 -(η + ū) 1/2 , (3.20)
and ρ(t) is given by equation (3.11).

4 Numerical solutions illustrating the results of Propositions 1, 2, 3

In this section, we construct numerical solutions to the nonlocal parabolic equation (2.2) complemented with the initial condition (3.7). In Subsection 4.1, we describe the setup of numerical simulations and the numerical methods employed. In Subsections 4.2-4.4, we present a sample of numerical solutions that illustrate the results established by Propositions 1, 2 and 3. In Subsection 4.5, we discuss the biological insights provided by these results.

Numerical methods and setup of numerical simulations

The method for constructing numerical solutions to the nonlocal parabolic equation (2.2) is based on a time splitting scheme between the conservative part and the reaction term. In particular, we approximate the diffusion term through a three-point finite difference explicit scheme and we use an implicit finite difference scheme for the reaction term [START_REF] Leveque | Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems[END_REF]. We select a uniform discretisation consisting of 500 points on the interval [-5, 5] as the computational domain of the independent variable x and we discretise the time interval [0, T ] with the uniform step ∆t = 0.001. On the other hand, we approximate the system of nonlinear ODEs (3.9) for f (t), µ(t) and ρ(t) by using an explicit Euler method with time-step ∆t = 0.0025. Numerical simulations are done in Matlab using the parameter values listed in Table 1, which have been estimated based on existing data from in vitro experiments on the phenotypic evolution of HL60 leukemic cells exposed to vincristine [START_REF] Pisco | Nondarwinian dynamics in therapy-induced cancer drug resistance[END_REF]. We refer the interested reader to the work by Lorenzi et al. [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF] for a detailed description of the model calibration. To reproduce the biological scenario in which the cell population has never been exposed to therapy, we define the terms ρ 0 , f 0 and µ 0 in equation (3.7) equal to the equilibrium values ρ ∞ ū=0 , f ∞ ū=0 and µ ∞ ū=0 given by equations (3.2) and (3.3), that is,

ρ 0 = 1 d γ -(β η) 1/2 , f 0 = η, and µ 0 = 0. (4.1)

Numerical solutions illustrating the results of Proposition 1

The plots presented in Figure 1 show a comparison between the numerical solution n(x, t) of equation ( 2.2) complemented with the initial condition (3.7) and the exact solution N (x, t) given by equation (3.8), with f (t), µ(t) and ρ(t) obtained by solving numerically the Cauchy problem (3.9). The graph of the drug delivery schedule

u(t) = 2 + sin(t/π) + cos( √ 3t/π) (4.2)
used to carry out numerical simulations is displayed in the upper row, left panel of Figure 1. We chose this example of a drug delivery schedule, with two incommensurate periodic components, to provide a more severe test of the agreement of the two computational approaches than a simpler drug delivery schedule would. The plots in the upper row of Figure 1 show that, in agreement with the results established by Proposition 1, there is a perfect match between the population size obtained by computing the integral of the numerical solution n(x, t) (upper row, central panel, solid red line) and the population size ρ(t) obtained by solving the Cauchy problem (3.9) (upper row, central panel, dashed blue line), as well as between the mean phenotypic state computed from the numerical solution n(x, t) (upper row, right panel, solid red line) and the mean phenotypic state µ(t) obtained by solving the Cauchy problem (3.9) (upper row, right panel, dashed blue line). Moreover, the results displayed in the lower row of Figure 1 indicate that there is also an excellent agreement between the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7) and the exact solution N (x, t) given by equation (3.8).

Numerical solutions illustrating the results of Proposition 2

The plots presented in Figure 2 show a comparison between the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7) and the exact solution N (x, t) given by equation (3.8), with ρ(t), f (t) and µ(t) computed through equations (3.11), (3.13) and (3.14), respectively. The graph of the drug delivery schedule u(t) used to carry out numerical simulations, which is of the form given by equation (3.12), is displayed in the upper row, left panel of Figure 2. There is an excellent agreement between the numerical results and the results established by Proposition 2.

Numerical solutions illustrating the results of Proposition 3

With the aim of illustrating the results of Proposition 3, in Figure 3 we compare the numerical solution n(x, t) of equation (2.2) with the exact solution N (x, t) given by equation (3.8), with ρ(t), f (t) and µ(t) given by equations (3.11), (3.18) and (3.19), respectively. Simulations have been carried using a drug delivery schedule of the form given by equation (3.17) with ū = 5.79, which corresponds to the LD30 dose (vid. upper row, left panel of Figure 3). The results shown in Figure 3 indicate that there is a perfect match between our numerical solutions and the analytical results established by Proposition 3. As one would expect based on the results of Theorem 1, after an initial transient: the population density converges to the steady state distribution given by equation (3.3); the total number of cells converges to the asymptotic value (3.2); the mean phenotypic state converges to the asymptotic value given by equation (3.3). 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1). 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1). ), respectively (central panel). The colour scale ranges from blue (low density) to yellow (high density). The plot in the right panel shows the comparison between the numerical solution n(x, t) (solid red line) and the exact solution N (x, t) (dashed blue line) at the final time t = T . Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Biological insights into cancer cell population evolution under cytotoxic therapy

In anticancer therapy, it is common to administer a given dose of cytotoxic agents at the beginning of a therapy cycle and then let the patient recover from toxicities. For this reason, we consider here a periodic, piecewiseconstant drug delivery schedule whereby a fixed drug dose is administered for 4 days followed by a rest period of 3 days over a time window of 12 weeks (i.e. 84 days). Such a drug delivery schedule captures the essentials of standard administration protocols for cytotoxic agents that are used clinically to treat cancer [START_REF] Skeel | Handbook of cancer chemotherapy[END_REF]. Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

The cytotoxic drug induces a population bottleneck

The numerical results shown in Figure 4 reveal that the cytotoxic drug induces a population bottleneck. In more detail, since the initial cell population is mainly composed of fast-proliferating cells in the phenotypic state x = 0 (i.e. cells which are sensitive to the cytotoxic drug) a drastic reduction of the population size occurs during the first therapy cycle (vid. central panel of Figure 4). This is followed by a short period of relatively constant population level, which lasts until cells in other phenotypic states closer to x = 1 (i.e. the phenotypic state corresponding to the highest level of cytotoxic-drug resistance) are selected and start proliferating, thus bringing about a sharp increase in the mean phenotypic state of the cell population (vid. right panel of Figure 4). As soon as the number of drug-tolerant cells is sufficiently large (i.e. the mean phenotypic state is sufficiently close to 1), the population becomes sufficiently resistant to the cytotoxic drug and the total number of cells regrows. From then on, successive therapy cycles produce only fluctuations of the cell population size and the mean phenotypic state about some stable values. Analogous considerations hold for the results presented in Figures 123. Moreover, in analogy with the results shown in Figures 123, the phenotypic distribution of the cell population remains unimodal throughout the entire time, with the mean phenotypic state being at the maximum point of the distribution (data not shown because qualitatively similar to those presented in the lower rows of Figures 123). Higher drug doses correspond to lower cell numbers but higher levels of drug resistance

The plots in Figure 5 depict the dynamics of the population size and the mean phenotypic state for the same simulation setup of Figure 4, but with different values of the fixed drug dose delivered during the administration phases. The curves in Figure 5 are qualitatively similar to those presented in Figure 4, but they also highlight how higher drug doses correspond to lower cell numbers and higher levels of drug resistance. In particular, the curves displayed in the left panel of Figure 5 show that, as one would expect, the cytotoxic drug has a detrimental effect on the population size. On the other hand, the curves displayed in the right panel of Figure 5 support the idea that different doses of the cytotoxic drug may lead the cell population to take different evolutionary trajectories, which culminate in the selection for phenotypic variants characterised by different levels of drug tolerance.

How the average population size varies with the evolutionary parameters of the model

The heatmaps in Figure 6 illustrate how the average size of the cancer cell population varies as a function of the rate of spontaneous epimutations β and of the selection gradient η (left panel), and as a function of the rate of death due to intrapopulation competition d and the maximum fitness γ (right panel). These results are for the drug delivery schedule shown in the left panel of Figure 4, but similar results have been obtained for different delivery schedules (data not shown). The average population size is a decreasing function of both β and η, which suggests that more frequent spontaneous epimutations and stronger selection pressures lead to cancer cell populations of a smaller size. Moreover, the average population size is an increasing function of the maximum fitness γ and a decreasing function of the average rate of cell death due to intrapopulation competition d. This is in line with the biological observation that cancer cell populations shrink in the presence 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

of harsher environments (i.e. lower values of γ) and more intense competition for limited resources (i.e. higher values of d).

The cancer cell population regains sensitivity to the cytotoxic drug during 'drug-holiday'

As illustrated by the results presented in Figure 7, when the drug treatment is interrupted, the mean phenotypic state decreases towards the weakly resistant and highly proliferating state x = 0 (vid. right panel of Figure 7) due to the fact that the number of more rapidly proliferating and drug-sensitive cells increases. Consequently, the population size regrows (vid. central panel of Figure 7). These results communicate the biological notion that while cells in a more resistant phenotypic state detain a competitive advantage against fast-proliferating cells in the presence of the cytotoxic drug, selection is reverted in favour of drug-sensitive cells with a higher proliferation rate during 'drug-holidays' [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF][START_REF] Thakur | Modelling vemurafenib resistance in melanoma reveals a strategy to forestall drug resistance[END_REF].

Adaptation to the cytotoxic drug is driven by the interplay between the drug selective pressure and selective forces that preexist the drug administration

The heatmaps in Figure 8 illustrate how the time average of the mean phenotypic state of cancer cells varies as a function of the rate of spontaneous epimutations β and of the selection gradient η. The results in the left panel are for the drug delivery schedule shown in Figure 4, whereas those in the right panel are for the drug delivery schedule of Figure 7. These results show that the time average of the mean phenotypic state is a decreasing function of the selection gradient η. Combining this fact together with the results presented in the right panel of Figure 5, one can reach the conclusion that cancer cell adaptation to the cytotoxic drug results Graph of the piecewise-constant drug delivery schedule u(t) used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

from the concerted action between selective forces that preexist the drug administration -the strength of which is measured by the parameter η -and the selective pressure exerted by the drug -which increases with the drug concentration. In fact, for the same value of the selection gradient, higher drug concentrations lead to values of the time average of the mean phenotypic state that are closer to the highly-resistant state x = 1. On the other hand, for the same drug concentration, higher values of the selection gradient lead to lower values of the time average of the mean phenotypic state, which correspond to lower levels of resistance and higher proliferation rates. Furthermore, the time average of the mean phenotypic state is a decreasing function of β. This is more evident in the presence of longer treatment breaks (vid. right panel of Figure 8) and is due to the fact that, in the absence of the drug, the mean phenotypic state of the cell population decreases more quickly towards the weakly resistant and highly proliferating state x = 0 for larger values of β.

Numerical optimal control

In the framework of our model, the problem of finding optimal dosage regimens that allows one to minimise the size of the cancer cell population can be tackled as an optimal control problem whereby the total number of cells ρ(t) and the drug dose u(t) represent, respectively, the state variable and the control variable. This is the subject of the present section. In Subsection 5.1, we formulate the optimal control problem, while in Subsection 5.2 we describe the method used to solve numerically such an optimal control problem and summarise the setup of numerical simulations. In Subsections 5.3-5.5, we present a sample of the results obtained and discuss their possible biological and clinical implications.

Optimal control problem

We consider both the biological scenario where one wishes to minimise the average number of cancer cells during the course of treatment (i.e. the mean value of the population size ρ(t) over the time interval [0, T ]) and the biological scenario where one's goal is to minimise the number of cancer cells at the end of the treatment (i.e.

Figure 8: Adaptation to the cytotoxic drug is driven by the interplay between the drug selective pressure and selective forces that preexist the drug administration. Time average of the mean phenotypic state µ(t) as a function of the rate of spontaneous epimutations β (in units of 10 -3 ) and of the selection gradient η, for the drug delivery schedules shown in Figure 4 (left panel) and Figure 7 (right panel). These numerical results have been obtained using a time window corresponding to 84 days with the values of the other model parameters as in Table 1 and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

the population size ρ(t) at the final time t = T ). Therefore, we define the cost functional either as

J = 1 T T 0 ρ(t) dt (5.1)
or as J = ρ(T ).

(5.2)

Several of the anticancer therapeutic protocols based on the administration of cytotoxic agents rely on dosage regimens which can be seen, in a first approximation, as piecewise-constant functions of time [START_REF] Skeel | Handbook of cancer chemotherapy[END_REF]. For this reason, we let the control u(t) be of the form given by equation (3.12), that is,

u(0) = u 0 , u(t) = u i 1 (τi,τi+1] (t) for i = 0, . . . , N -1, (5.3)
where

τ 0 = 0 < τ 1 < . . . < τ N -1 < τ N = T and u = (u 0 , . . . , u N -1 ) ∈ R N ≥0 , |u| ∞ < u LD100 , (5.4)
with u LD100 ∈ R >0 being the smallest value of ū such that under the constant drug therapy u(t) = ū we have ρ(t) → 0 as t → ∞, as discussed in Section 3 (vid. Remark 2). We assume the values of the parameters τ 0 , . . . , τ N to be given and define the space of admissible controls U as the following finite-dimensional space

U = u ∈ R N ≥0 : |u| ∞ < u LD100 . (5.5)
In this setting, we consider the optimal control problem inf u∈U J,

for J defined either by equation (5.1) or by equation (5.2), subject to the dynamic constraints given by the state equation

ρ (t) = R R(x, ρ(t), u(t)) n(x, t) dx, (5.7)
with n(x, t) being the solution of the nonlocal parabolic equation (2.2) complemented with the initial condition (3.7), and the control constraints

k+M k j=k u j ≤ C k with C k ∈ R ≥0 , M k ∈ [0, N -1] ⊂ N 0 and k ∈ [0, N -1 -M k ] ⊂ N 0 . (5.8)
The constraints given by equation (5.8) translate into mathematical terms the idea that, due to the documented side effects of cytotoxic agents, one may want the drug dose administered during given time intervals to remain below certain maximum levels. The parameter M k models the duration of the therapy cycle labelled by the index k. Moreover, the parameter C k represents the maximum tolerated dose (MTD) for the therapy cycle k, which corresponds to the highest total dose that can be administered over the therapy cycle without causing unacceptable side effects. If M 0 = N -1, then the only admissible value of k is 0 and, after having renamed C 0 to C, one can rewrite the control constraint given by equation (5.8) as

|u| 1 ≤ C, (5.9)
with C being the maximum administrable dose during the course of the whole treatment (i.e. over the whole time interval [0, T ]).

Numerical methods and setup of numerical simulations

In order to reduce the computational cost of solving the optimal control problem under study, we exploit the analytical results established by Proposition 1, which allow one to replace the dynamic constraints given by the IDE (5.7) with the equivalent constraints given by the system of ODEs (3.9). Hence, we solve numerically the optimal control problem given by equation (5.6), for J defined either by equation (5.1) or by equation (5.2), subject to the constraints given by equations (3.9) and (5.8). Numerical solutions are constructed using the interior point method [START_REF] Nocedal | Numerical optimization[END_REF]. Numerical simulations are done in Matlab using the parameter values given in Table 1 and the initial condition (3.7), with the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1). Moreover, in equation (5.4) we define τ i = iτ for all i = 0, . . . , N -1 with the value of τ corresponding to 1 day.

(5.10)

The parameter choice given by equation (5.10) translates into mathematical terms the idea that, due to constraints associated with clinical feasibility, one cannot vary the delivered drug dose in the span of 1 day. However, we remark that our mathematical framework and the related numerical methods can accommodate any values of the parameter τ . As for the constraints given by equation (5.8), we consider different values of the parameters M k (i.e. the duration of the therapy cycle k) and C k (i.e. the MTD for the therapy cycle k) which correspond to different biological and clinical scenarios, and we define the values of C k in terms of the LDα.

Optimal dosing regimen for 12 therapy cycles of 1 week

We focus on a treatment protocol which consists of 12 therapy cycles of 1 week with a constraint on the total amount of cytotoxic drug that can be administered per cycle, which we assume to be the same for all the therapy cycles. For this reason, we choose a time window corresponding to 84 days (i.e. 12 weeks) and, under the assumptions given by equation (5.10), we specify the constraint given by equation (5.8) as follows

7 h+7
j=7 h u j ≤ C with C ∈ R >0 and h = 0, 1, . . . , 10, 11.

(5.11)

We tested different values of the MTD per therapy cycle C ranging from the LD30 to the LD80. For J defined by equation (5.1), the results obtained with values of C between the LD40 and the LD80 are qualitatively similar. Therefore, we report here on the outcomes of numerical simulations for C corresponding to the LD30 and the LD80. Moreover, for J defined by equation (5.2), the results obtained are qualitatively similar for all values of C. Therefore, we report here on the outcomes of numerical simulations for C corresponding to the LD80 only.

The numerical solutions presented in the upper row, left panel of Figure 9 show that, when the goal is to minimise the average number of cancer cells during the course of treatment and the MTD per therapy cycle corresponds to the LD30, the optimal dosing regimen is as follows: cycles 1-5: delivery of the MTD on the first day of the cycle followed by 6 days of rest; cycles 6-9: delivery of a higher dose on the first day of the cycle, which is monotonically decreasing from one cycle to the other, followed by a variable but consistently low drug dose distributed over a number of days -the number of days monotonically decreases from one cycle to the other; cycle 10: delivery of a higher dose on the first and the last day of the cycle interspersed with the administration of a variably lower drug dose; cycle 11: continuous delivery of a monotonically increasing, but consistently low, drug dose; cycle 12: delivery of a relatively low and variable drug dose during the first 3 days of the cycle followed by 4 days of rest. The amount of cytotoxic drug administered per therapy cycle is such that the constraint given by equation (5.11) is saturated for all values of h. Under this optimal dosing schedule, the number of cancer cells initially decreases (vid. upper row, central panel of Figure 9) but then grows again due to the increase of the mean level of drug resistance (vid. upper row, right panel of Figure 9). Ultimately, both the population size and the mean phenotypic state appear to oscillate about some stable values.

A clear change in the optimal dosing regimen occurs when higher values of the MTD per therapy cycle are allowed. In particular, the numerical results presented in the lower row, left panel of Figure 9 show that, when the MTD per therapy cycle is the LD80, the optimal drug delivery schedule is as follows: cycle 1: delivery of the MTD on the first day of the cycle followed by 6 days of rest; cycle 2: delivery of a higher dose on the first day of the cycle followed by a variable, but consistently low, drug dose over days 2-6 and then one day of rest; cycles 3 and 4: delivery of a higher dose on the first day of the cycle, which is monotonically decreasing from one cycle to the other, followed by a variable but consistently low drug dose administered throughout the rest of the cycle; cycles 5-10: continuous delivery of a low and almost constant drug dose; cycle 11: continuous delivery of a variable, but consistently low, drug dose; cycle 12: delivery of a low and variable drug dose during the first 5 days of the cycle followed by 2 days of rest. The amount of cytotoxic drug administered per cycle is such that the constraint given by equation (5.11) is saturated for all values of h. Under this optimal dosing schedule, the dynamics of the number of cancer cells and of the mean phenotypic state displayed in the lower row of Figure 9 are qualitatively similar to those shown in the upper row. However, in this case the oscillations are much less pronounced. Moreover, since higher drug doses are allowed in each therapy cycle, the mean population size and the time average of the mean phenotypic state are, respectively, lower and higher than those of the results presented in the upper row (cf. the values highlighted by the dotted black lines in the central and right panels of Figure 9).

On the other hand, in the case where the goal is to minimise the number of viable cancer cells at the end of the treatment, the numerical solutions presented in the left panel of Figure 10 suggest that the optimal dosing regimen consists in administering the MTD on the final day of the last therapy cycle. Therefore, the Figure 9: Optimal dosing regimen for 12 therapy cycles of 1 week with the goal being to minimise the average total number of cells during the course of treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.1), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1) and the additional assumptions given by equation (5.10). Focussing on a treatment protocol which consists of 12 therapy cycles of 1 week with a constraint on the total amount of drug that can be administered per cycle, which is the same for all the therapy cycles, we have specified the constraint given by equation ( 5 constraint given by equation (5.11) is saturated only for h = 11. Under such an optimal dosing schedule, the number of cancer cells is constant for all the time but on the last day of treatment (i.e. on day 84), when it undergoes a drastic decay -vid. central panel of Figure 10. This is due to the fact that the initial number of cells corresponds to the equilibrium population size without drug [i.e. the value ρ ∞ ū=0 given by equation (3.2)]. If lower initial cell numbers are considered, then the population size grows until it saturates at the value ρ ∞ ū=0 from day 0 to day 83 included and decays on day 84, when the cytotoxic drug is administered (vid. Figure S1 in Appendix E). Moreover, the mean phenotypic state remains constant before undergoing a sharp increase on day 84, due to the occurrence of directional selection for resistant phenotypic variants, which is induced by the delivery of the cytotoxic drug (vid. right panel of Figure 10). Figure 10: Optimal dosing regimen for 12 therapy cycles of 1 week with the goal being to minimise the total number of cells at the end of the treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation ( 5.2), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms ρ 0 , f 0 and µ 0 given in equation ( 4.1) and the additional assumptions given by equation (5.10). Focussing on a treatment protocol which consists of 12 therapy cycles of 1 week with a constraint on the total amount of drug that can be administered per cycle, which is the same for all the therapy cycles, we have specified the constraint given by equation (5.8) as in equation (5.11) with the value of C corresponding to the LD80. The coloured rectangles in the left panel highlight the different therapy cycles.

Optimal dosing regimen for 12 therapy cycles of 1 week with 3 days of rest per week

We consider a treatment protocol which consists of 12 therapy cycles of 1 week whereby 4 days of therapy are followed by 3 days of rest. Furthermore, we assume that there is a maximum amount of drug that can be administered over the 4 days of therapy, which is the same for all the therapy cycles. Hence, we choose a time window corresponding to 84 days (i.e. 12 weeks) and, under the assumptions given by equation (5.10), we specify the constraint given by equation (5.8) as follows u j = 0 with C ∈ R >0 and h = 0, 1, . . . , 10, 11.

(5.12)

We tested different values of the MTD per therapy cycle C ranging from the LD30 to the LD80. For J defined by equation (5.1), the results obtained with values of C between the LD40 and the LD80 are qualitatively similar. Therefore, we report here on the outcomes of numerical simulations for C corresponding to the LD30 and the LD80. Moreover, for J defined by equation (5.2), the results obtained are qualitatively similar for all values of C. Therefore, we report here on the outcomes of numerical simulations for C corresponding to the LD80 only.

The numerical solutions presented in the upper row, left panel of Figure 11 show that, when the goal is to minimise the average number of cancer cells during the course of treatment and the MTD per therapy cycle corresponds to the LD30, the optimal dosing regimen is as follows: cycles 1-5: delivery of the MTD on the first day of the cycle followed by 6 days of rest; cycles 6-9: delivery of a higher dose on the first day of the cycle, which is monotonically decreasing from one cycle to the other, followed by a variable but consistently low drug dose distributed over days 2-4; cycles 10 and 11: delivery of a higher dose on day 4 and lower variable drug doses on the other days; cycle 12: delivery of a relatively low and decreasing drug dose during days 1-4. The amount of cytotoxic drug administered per cycle is such that the constraint given by equation (5.12) is saturated for all values of h. Under this optimal dosing schedule, the number of cancer cells initially decreases (vid. upper row, central panel of Figure 11) but then grows again due to the increase of the mean level of drug resistance (vid. upper row, right panel of Figure 11). Ultimately, both the population size and the mean phenotypic state appear to oscillate about some stable values.

A clear change in the structure of the optimal dosing regimen occurs when higher values of the MTD per therapy cycle are considered. In particular, the numerical results presented in the lower row, left panel of Figure 11 show that, when the MTD per therapy cycle is the LD80, the optimal drug delivery schedule is as follows: cycle 1: delivery of the MTD on the first day of the cycle; cycles 2-12: delivery of a moderately high dose on day 1 and day 4 interspersed with the continuous administration of a relatively low drug dose. The amount of cytotoxic drug administered per cycle is such that the constraint given by equation (5.12) is saturated for all values of h. Under this optimal dosing schedule, the sample dynamics of the number of cancer cells and of the mean phenotypic state displayed in the lower row of Figure 11 are qualitatively similar to those shown in the upper row. However, in this case the oscillations in the mean phenotypic state are much less pronounced. Moreover, since higher drug doses are allowed in each therapy cycle, the mean population size and the time average of the mean phenotypic state are, respectively, lower and higher than those of the results presented in the upper row (cf. the values highlighted by the dotted black lines in the central and right panels of Figure 11).

On the other hand, in the case where the goal is to minimise the number of cancer cells at the end of the treatment, the numerical solutions presented in the left panel of Figure 12 show that the optimal dosing regimen consists in administering the MTD on the fourth day of the last therapy cycle (i.e. the last day on which the cytotoxic drug can be delivered). Therefore, the constraint given by equation (5.11) is saturated only for h = 11. Under such an optimal dosing schedule, the number of cancer cells is constant for most of the time, decreases drastically on the fourth day of the last therapy cycle (i.e. day 80) and then starts increasing slowly during the last phase of imposed rest (i.e. days 81-84) -vid. central panel of Figure 12. This is due to the fact that the initial number of cells corresponds to the equilibrium population size without drug [i.e. the value ρ ∞ ū=0 given by equation (3.2)]. If lower initial cell numbers are considered, then the population size grows until it saturates at the value ρ ∞ ū=0 from day 0 to day 79 included and decays on day 80, when the cytotoxic drug is administered (vid. Figure S2 in Appendix E). Moreover, the mean phenotypic state remains constant before undergoing a sharp increase on day 80, due to the occurrence of directional selection for drug resistant phenotypic variants, and then starts decreasing slowly during days 81-84 (vid. right panel of Figure 12).

Comparing the administration of high drug doses separated by drug-free periods with the continuous administration of a relatively low drug dose

In this subsection, we compare the therapeutic performance of the periodic administrations of high drug doses separated by drug-free periods with the therapeutic performance of the continuous administration of a relatively Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.1), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1) and the additional assumptions given by equation (5.10). The simulations refer to a treatment protocol which consists of 12 therapy cycles of 1 week whereby 4 days of therapy must be followed by 3 days of rest, and there is a maximum amount of drug that can be administered over the 4 days of therapy, which is the same for all the therapy cycles. Therefore, we have specified the constraint given by equation (5.8) as in equation (5.12) with the value of C corresponding either to the LD30 (upper row) or to the LD80 (lower row). The coloured rectangles in the left panel highlight the different therapy cycles. The dotted black lines in the central and right panels highlight the values of the time averages of the total number of cells and of the mean phenotypic state.

low drug dose, in relation to the performance of the optimal dosing regimens presented in the lower row of Figure 9 and in Figure 10 -i.e. we focus on a therapy protocol which consists of 12 therapy cycles of 1 week with a given MTD per therapy cycle, which is the same for all the cycles and corresponds to the LD80. The therapeutic performance is assessed in terms of either the average cell population size during the course of treatment or the size of the cancer cell population at the end of the treatment. The aim is to verify which between the administration of high drug doses separated by drug-free periods and the continuous administration Figure 12: Optimal dosing regimen for 12 therapy cycles of 1 week with 3 days of rest per week with the goal being to minimise the total number of cells at the end of the treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.2), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8).

Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1) and the additional assumptions given by equation (5.10). The simulations refer to a treatment protocol which consists of 12 therapy cycles of 1 week whereby 4 days of therapy must be followed by 3 days of rest, and there is a maximum amount of drug that can be administered over the 4 days of therapy, which is the same for all the therapy cycles. Therefore, we have specified the constraint given by equation (5.8) as in equation (5.12) with the value of C corresponding to the LD80. The coloured rectangles in the left panel highlight the different therapy cycles.

of a relatively low drug dose performs more closely to the optimal dosing regimen. For the administration of high drug doses separated by drug-free periods, we consider 12 therapy cycles of 1 week whereby the MTD is administered either on day 1, day 4 or day 7 of every cycle (vid. blue, green or yellow line in the left panel of Figure 13). As for the continuous administration of a relatively low drug dose, we consider the constant delivery for 84 days (i.e. 12 weeks) of an amount of drug such that the total dose administered over the whole treatment is the same as that of the other drug delivery schedules (vid. cyan line in the left panel of Figure 13). The numerical results obtained in the case where one wishes to minimise the average population size during the course of treatment, which are summarised by the plot in the central panel of Figure 13, indicate that the continuous administration of a relatively low drug dose performs closely to the optimal dosing regimen, and better than the other delivery schedules. On the other hand, as illustrated by the plot in the right panel of Figure 13, in the case where the goal is to minimise the size of the cell population at the end of the treatment, the outcomes of our model support the idea that the best performance is achieved by the drug delivery schedule whereby the MTD is delivered on the the 7th day of every cycle. Notice that in this case there is a relatively wide spread between the performance of the optimal dosing regimen and that of the other drug delivery schedules under consideration. ) obtained under the optimal dosing regimen shown in the lower row, left panel of Figure 9 [i.e. u opt (t)] and the drug delivery schedules u 1 (t), u 2 (t), u 3 (t) and u 4 (t). Right panel. Comparison between the final cell population size (in units of 10 8 ) obtained under the optimal dosing regimen shown in the left panel of Figure 10 [i.e. u opt (t)] and the drug delivery schedules u 1 (t), u 2 (t), u 3 (t) and u 4 (t).

Conclusions

In this paper, we have constructed exact solutions and addressed numerical optimal control for a nonlocal parabolic equation modelling evolutionary dynamics in a population of cancer cells exposed to the action of a cytotoxic drug. The results that we have obtained give answers to the open questions (i)-(iii) posed in the introduction of this paper. In summary:

(i) Our results indicate that the cytotoxic drug can cause a drastic reduction in the size of the cancer cell population during the first therapy cycles, when most of the cells are in a fast-proliferating and drugsensitive phenotypic state. However, due to the selective sweep towards drug-tolerant phenotypic variants induced by the cytotoxic drug, the cell population progressively acquires drug resistance. As a consequence, the population size regrows after some cycles of therapy. From then on, successive administrations of the cytotoxic drug produce only fluctuations of the cell population size about some stable values.

(ii) In the case where one wishes to minimise the number of viable cancer cells at the end of the treatment, our results show that a therapeutic performance close to that of the optimal dosing regimen can be achieved by the periodic administration of high drug doses separated by drug-free periods.

(iii) When the desired goal is to minimise the average size of the cancer cell population during the course of treatment, our results suggest that the continuous administration of a relatively low drug dose can give similar performance to the optimal dosing regimen.

Our theoretical work illustrates how cancer cell adaptation to cytotoxic therapy results from a concerted action between the selective pressure exerted by the drug in use and preexisting selective forces which are in place within the cell population [START_REF] Moreno | Is cell competition relevant to cancer?[END_REF][START_REF] Pienta | Ecological therapy for cancer: defining tumors using an ecosystem paradigm suggests new opportunities for novel cancer treatments[END_REF][START_REF] Tamori | Cell competition and its implications for development and cancer[END_REF][START_REF] Wagstaff | Competitive cell interactions in cancer: a cellular tug of war[END_REF]. Furthermore, in agreement with previous experimental results [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF][START_REF] Sun | Reversible and adaptive resistance to braf (v600e) inhibition in melanoma[END_REF], our model suggests that cancer cell populations may regain sensitivity to cytotoxic drugs after a 'drug holiday'. The results obtained also show that the level of drug resistance acquired by the cell population can vary with the administered drug dose, which suggests that different doses of the same cytotoxic agent can trigger the selection for phenotypic variants characterised by different levels of drug tolerance. This supports the notion that the response of cancer cells to cytotoxic agents is more complex than a simple binary outcome -i.e. extinction of sensitive cells and selection for highly resistant cells -and recapitulates the results from previous theoretical and experimental studies [START_REF] Bozic | Dynamics of targeted cancer therapy[END_REF][START_REF] Burrell | Tumour heterogeneity and the evolution of polyclonal drug resistance[END_REF][START_REF] De Bruin | Intra-tumor heterogeneity: lessons from microbial evolution and clinical implications[END_REF][START_REF] Souza | Chemotherapy dosing schedule influences drug resistance development in ovarian cancer[END_REF][START_REF] Enriquez-Navas | Exploiting evolutionary principles to prolong tumor control in preclinical models of breast cancer[END_REF][START_REF] Gatenby | Adaptive therapy[END_REF][START_REF] Gatenby | A change of strategy in the war on cancer[END_REF][START_REF] Merlo | Cancer as an evolutionary and ecological process[END_REF][START_REF] Salgia | The genetic/non-genetic duality of drug 'resistance'in cancer[END_REF].

In current clinical practice, it is common to give the MTD of cytotoxic agents at the beginning of a therapy cycle and then let the patient recover from toxicities. In the last decade, experimental and theoretical efforts have been made to determine the best scheduling of cytotoxic agents in anticancer therapy; unfortunately, there is still no clear answer about optimal administration protocols [START_REF] Savage | Why does cytotoxic chemotherapy cure only some cancers?[END_REF]. Our theoretical work demonstrates that higher doses of cytotoxic drugs reduce the size of cancer cell populations at the cost of promoting the selection for more resistant phenotypic variants. Moreover, our results provide a conceptual basis for the benefits of metronomic therapy [START_REF] André | Metronomics: towards personalized chemotherapy?[END_REF][START_REF] Benzekry | Metronomic reloaded: theoretical models bringing chemotherapy into the era of precision medicine[END_REF][START_REF] Hanahan | Less is more, regularly: metronomic dosing of cytotoxic drugs can target tumor angiogenesis in mice[END_REF][START_REF] Klement | Continuous low-dose therapy with vinblastine and vegf receptor-2 antibody induces sustained tumor regression without overt toxicity[END_REF][START_REF] Munzone | Clinical overview of metronomic chemotherapy in breast cancer[END_REF][START_REF] Scharovsky | Metronomic chemotherapy: changing the paradigm that more is better[END_REF]] -which relies on continuous treatment with lower doses of cytotoxic agents -as it performs more closely to the optimal dosing regimen in the case where the goal is to minimise the average number of viable cancer cells during the course of treatment. The parameter values that we have used to carry out numerical simulations are derived from leukemia datasets. However, given the robustness and structural stability of our results, we expect these conclusions to hold for different types of cancer.

We conclude with an outlook on possible extensions of the present work. A natural way to extend our work would be to introduce an advection term to take into account the effects of stress-induced epimutations, which are triggered by the selective pressure that cytotoxic drugs exert on cancer cells [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF]. The present study could also be extended by including a spatial structure. For instance, one could embed the cancer cells in the geometry of an avascular (or a vascular) tumour and introduce an equation for the evolution of the local concentration of the cytotoxic drug in the interior of the tumour [START_REF] Lorenzi | The role of spatial variations of abiotic factors in mediating intratumour phenotypic heterogeneity[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]. As for the optimal control problem, in such a spatial setting one could select the total number of viable cells within the tumour as the state variable and the drug concentration at the boundary of the tumour (or in the blood vessels) as the control variable. Finally, since continuum models like the one that we have considered can be derived from stochastic individual-based models [START_REF] Champagnat | The canonical equation of adaptive dynamics: a mathematical view[END_REF][START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF][START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to population-level consequences[END_REF][START_REF] Stace | A phenotype-structured individual-based model for the evolution of cancer cell populations under chemotherapy[END_REF], which provide a more accurate modelling framework in the regime of low cell numbers, another track to follow might be to investigate whether the drug delivery schedules obtained through the numerical optimal control of the continuum model would remain optimal also for an equivalent individual-based model.

with

ρ(t) = R N (x, t) dx and µ(t) = 1 ρ(t) R x N (x, t) dx. Since log (N ) = log (ρ(t)) + 1 4 log f (t) 4β - f (t) 4β 1/2 (x -µ(t)) 2 + constant, it follows that 1 N ∂N ∂t = ρ (t) ρ(t) + 1 4
f (t) f (t) - 1 4
f (t) (βf (t)) 1/2 (x -µ(t)) 2 + f (t) β 1/2 (x -µ(t))µ (t), 1 N ∂N ∂x = - f (t) β 1/2 (x -µ(t)), 1 N ∂ 2 N ∂x 2 = - f (t) β 1/2 + f (t) β (x -µ(t)) 2 .
Hence, substituting the ansatz N (x, t) into equation (A.2) one finds that

ρ (t) ρ(t) + 1 4
f (t) f (t) - 1 4
f (t) (βf (t)) 1/2 [x -µ(t)] 2 + f (t) β 1/2 (x -µ(t))µ (t) = β - f (t) β 1/2 + f (t) β (x -µ(t)) 2 + Γ(u(t)) -H(u(t)) (x -Θ(u(t))) 2 -dρ(t).
Equating the coefficients of the zero-order, first-order and second-order terms in x produces a system of ODEs.

The second-order terms in x yield a differential equation for f alone, namely

f (t) = 4 H(u(t)) (βf (t)) 1/2 -4β 1/2 (f (t)) 3/2 . (A.4)
Equating coefficients of the first-order terms in x and eliminating f from the resulting equation produces

µ (t) = 2 β f (t) 1/2 H(u(t)) (Θ(u(t)) -µ(t)) . (A.5)
Finally, equating the zero-order terms in x and eliminating f and µ we find

ρ (t) ρ(t) = Q(t) -dρ(t), (A.6)
with Q(t) being given by equation (3.10), that is,

Q(t) = Γ(u(t)) -H(u(t)) (Θ(u(t)) -µ(t)) 2 -H(u(t)) β f (t) 1/2 .
If we impose the initial conditions f (0) = f 0 , µ(0) = µ 0 and ρ(0) = ρ 0 on the system of equations (A.4)-(A.6) to enforce the initial condition (3.7) for n(x, 0), we arrive at the Cauchy problem (3.9) for the functions f (t), µ(t) and ρ(t). Finally, using the substitution

ρ(t) = 1 d r(t) , ρ (t) ρ(t) = - r (t) r(t) (A.7)
one finds that r(t) satisfies the Cauchy problem r

(t) = 1 -Q(t)r(t), r(0) = [dρ(0)] -1 . Thus r(t) = exp - t 0 Q(s) ds 1 dρ(0) + t 0 exp z 0 Q(s)

B Proof of Proposition 2

Building upon the results of Proposition 1, in order to prove the claim of Proposition 2 it suffices to prove that the functions f (t) and µ(t) are given by the recursive schemes of equations (3.13)- (3.16). We note that Proposition 1 shows that Gaussian profiles for n(x, t) are preserved under time evolution. We evolve the system over the first interval of constant u(t), that is from t = τ 0 = 0 to t = τ 1 , and use the values of f (τ 1 ), µ(τ 1 ) and ρ(τ 1 ) as initial data for a new Cauchy problem for t ∈ [τ 1 , τ 2 ], and so on. Therefore, to complete the proof, we need only establish the solution under constant u(t) for one finite time interval (τ i , τ i+1 ]. For brevity, in this time interval we use compact versions of the notation given in equation (A.1):

Γ(u(t)) ≡ Γ i = γ - η u i η + u i , H(u(t)) ≡ H i = η + u i , Θ(u(t)) ≡ Θ i = u i η + u i .
It is also convenient to work with a time variable with initial value 0. Hence, we write

f i (t) = f (t + τ i ), t ∈ (0, τi ], τi = τ i+1 -τ i ,
and begin by solving the Cauchy problem for f i (t) for t ∈ (0, τi ]. In this interval we also write

g i (t) = (βf i (t)) 1/2 , (B.1)
and note that g i (0) = (βf i (0)) 1/2 = (βf (τ i )) 1/2 . We find that g i (t) evolves subject to the differential equation

g i (t) = 2(β H i -g 2 i (t)). (B.2)
This equation is of Riccati type and can be converted to a linear ODE of second order by an appropriate transformation [START_REF] Euler | De integratione aequationum differentialium[END_REF]. Writing

g i (t) = 1 2 d dt log(h i (t)) = h i (t) 2h i (t) (B.3) we find that h i (t) = 4βH i h i (t) (B.4)
and so

h i (t) = A cosh 2β 1/2 H 1/2 i t + B sinh 2β 1/2 H 1/2 i t , t ∈ (0, τi ], A and B constant. (B.5)
To ensure that h i (0) > 0 in the vicinity of t = 0 we need A > 0, while to ensure that lim t→0 + g i (t) > 0 we need h i (0) > 0, and so B > 0. We write C = B/A, so C ∈ (0, ∞). From equation (B.5) we find that

f i (t) = 1 4β h i (t) h i (t) 2 = H i   tanh 2β 1/2 H 1/2 i t + C 1 + C tanh 2β 1/2 H 1/2 i t   2 = H i   1 + C -1 tanh 2β 1/2 H 1/2 i t C -1 + tanh 2β 1/2 H 1/2 i t   2 (B.6)
and letting t = 0 shows that

C = f i (0) H i 1/2 . (B.7) If C = 1, that is, if f i (0) = H i , then f i (t) = H i for all t ∈ [0, τi ]. (B.8)
Therefore, we need to consider further only two cases, namely C < 1 and C > 1 and for each of these we can make use of the identity

tanh(z 1 + z 2 ) = tanh(z 1 ) + tanh(z 2 ) 1 + tanh(z 1 ) tanh(z 2 ) . (B.9) Case C < 1, that is, f i (0) < H i . If C < 1 we can write C = tanh arctanh (f i (0)/H i ) 1/2
and so we arrive at

f i (t) = H i tanh 2 2β 1/2 H 1/2 i t + arctanh f i (0) H i 1/2 , for t ∈ [0, τi ]. (B.10) Case C > 1, that is, f i (0) > H i . If C > 1 we can write C -1 = tanh arctanh (H i /f i (0)) 1/2
and we find that

f i (t) = H i coth 2 2β 1/2 H 1/2 i t + arctanh H i f i (0) 1/2 , for t ∈ [0, τi ]. (B.11)
Taken together, the results given by equations (B.8), (B.10) and (B.11) allow one to obtain the recursive scheme given by equation (3.13).

To prove the results given by equation (3.14), with τi = τ i+1 -τ i as above, we write µ i (t) = µ(t + τ i ), for t ∈ (0, τi ], so that we have to solve a Cauchy problem for µ i (t) for t ∈ (0, τi ], with µ i (0) = µ(τ i ). Replacing f (t) in the system of equations (3.9) by f i (t), we find that the function µ i (t) evolves under the differential equation In particular, choosing ξ = 0 we find that, where K is a positive constant, .

µ i (t) + 2 β f i (t) 1/2 H i µ i (t) = 2 β f i (t)
H i t 0 2 β f i (s)
We note that, since sinh (arctanh(z)) = tanh (arctanh(z)) cosh (arctanh(z)) = tanh (arctanh(z))

sech 2 (arctanh(z)) = z √ 1 -z 2 ,
we have

ζ i = (f i (0)/H i ) 1/2 (1 -f i (0)/H i ) 1/2 = f i (0) H i -f i (0) 1/2 .
We can obtain from equation (B.15) an explicit solution for µ i (t) when f i (0) < H i : In particular, choosing ξ = 0 we find that, where L is a positive constant, Taken together, these results allow one to obtain a recursive scheme for Q(t) and t 0 Q(s)ds. 2), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms f 0 and µ 0 given in equation (4.1), ρ 0 = ρ ∞ ū=0 /4 with ρ ∞ given by equation (3.2), and the additional assumptions given by equation (5.10). The simulations refer to a treatment protocol which consists of 12 therapy cycles of 1 week whereby 4 days of therapy must be followed by 3 days of rest, and there is a maximum amount of drug that can be administered over the 4 days of therapy, which is the same for all the therapy cycles. Therefore, we have specified the constraint given by equation (5.8) as in equation (5.12) with the value of C corresponding to the LD80. The coloured rectangles in the left panel highlight the different therapy cycles.

µ i (t) = Θ i + ζ i (µ i (0) -Θ i ) cosech 2β
H i t 0 2 β f i (s)

Figure 1 :

 1 Figure 1: Numerical solutions illustrating the results of Proposition 1. Upper row. Graph of the drug delivery schedule u(t) given by equation (4.2) that has been used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel). The solid red lines correspond to the quantities computed from the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7), while the dashed blue lines highlight ρ(t) and µ(t) computed by solving numerically the Cauchy problem (3.9). Lower row. Dynamics of the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7) (left panel) and of the exact solution N (x, t) given by equation (3.8) with f (t), µ(t) and ρ(t) obtained by solving numerically the Cauchy problem (3.9) (central panel). The colour scale ranges from blue (low density) to yellow (high density). The plot in the right panel shows the comparison between the numerical solution n(x, t) (solid red line) and the exact solution N (x, t) (dashed blue line) at the final time t = T . Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Figure 2 :

 2 Figure 2: Numerical solutions illustrating the results of Proposition 2. Upper row. Graph of the piecewise-constant drug delivery schedule u(t) that has been used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel). The solid red lines correspond to the quantities computed from the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7), while the dashed blue lines highlight ρ(t) and µ(t) computed through equations (3.11), (3.13) and (3.14). Lower row. Dynamics of the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7) (left panel) and of the exact solution N (x, t) given by equation (3.8) with f (t), µ(t) and ρ(t) computed through equations (3.11), (3.13) and (3.14) (central panel).The colour scale ranges from blue (low density) to yellow (high density). The plot in the right panel shows the comparison between the numerical solution n(x, t) (solid red line) and the exact solution N (x, t) (dashed blue line) at the final time t = T . Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Figure 3 :

 3 Figure 3: Numerical solutions illustrating the results of Proposition 3. Upper row. Graph of the constant drug delivery schedule u(t) that has been used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel). The solid red lines correspond to the quantities computed from the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7), while the dashed blue lines highlight ρ(t) and µ(t) computed through equations (3.11), (3.18) and (3.19). Lower row. Dynamics of the numerical solution n(x, t) of equation (2.2) complemented with the initial condition (3.7) (left panel) and of the exact solution N (x, t) given by equation (3.8) with ρ(t), f (t) and µ(t) given by equations (3.11), (3.18) and (3.19), respectively (central panel). The colour scale ranges from blue (low density) to yellow (high density). The plot in the right panel shows the comparison between the numerical solution n(x, t) (solid red line) and the exact solution N (x, t) (dashed blue line) at the final time t = T . Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Figure 4 :

 4 Figure 4: The cytotoxic drug induces a population bottleneck. Graph of the piecewise-constant drug delivery schedule u(t) that has been used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel).Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Figure 5 :

 5 Figure 5: Higher drug doses correspond to lower cell numbers but higher levels of drug resistance. Dynamics of the total number of cells ρ(t) in units of 10 8 (left panel) and of the mean phenotypic state µ(t) (right panel) for the same simulation setup of Figure 4 but with different values of the fixed drug dose delivered during the administration phases -i.e. LD10 (green curves), LD20 (blue curves), LD25 (cyan curves) and LD80 (yellow curves).

Figure 6 :

 6 Figure 6: How the average population size varies with the evolutionary parameters of the model. Mean value of the total number of cells ρ(t) in units of 10 8 as a function of the rate of spontaneous epimutations β (in units of 10 -3 ) and of the selection gradient η (left panel), and as a function of the rate of death due to intrapopulation competition d (in units of 10 -8 ) and of the maximum fitness γ (right panel). These numerical results are for the drug delivery schedule shown in the left panel of Figure 4. They have been obtained using a time window corresponding to 84 days with the values of the other model parameters as in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

Figure 7 :

 7 Figure 7: The cancer cell population regains sensitivity to the cytotoxic drug during 'drug-holiday'.Graph of the piecewise-constant drug delivery schedule u(t) used to carry out numerical simulations (left panel); dynamics of the total number of cells ρ(t) in units of 10 8 (central panel); dynamics of the mean phenotypic state µ(t) (right panel). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1and the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1).

  Figure9: Optimal dosing regimen for 12 therapy cycles of 1 week with the goal being to minimise the average total number of cells during the course of treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.1), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1, the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1) and the additional assumptions given by equation (5.10). Focussing on a treatment protocol which consists of 12 therapy cycles of 1 week with a constraint on the total amount of drug that can be administered per cycle, which is the same for all the therapy cycles, we have specified the constraint given by equation (5.8) as in equation(5.11) with the value of C corresponding either to the LD30 (upper row) or to the LD80 (lower row). The coloured rectangles in the left panel highlight the different therapy cycles. The dotted black lines in the central and right panels highlight the values of the time averages of the total number of cells and of the mean phenotypic state.

Figure 11 :

 11 Figure 11: Optimal dosing regimen for 12 therapy cycles of 1 week with 3 days of rest per week with the goal being to minimise the average total number of cells during the course of treatment.Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.1), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1, the values of the terms ρ 0 , f 0 and µ 0 given in equation (4.1) and the additional assumptions given by equation (5.10). The simulations refer to a treatment protocol which consists of 12 therapy cycles of 1 week whereby 4 days of therapy must be followed by 3 days of rest, and there is a maximum amount of drug that can be administered over the 4 days of therapy, which is the same for all the therapy cycles. Therefore, we have specified the constraint given by equation (5.8) as in equation (5.12) with the value of C corresponding either to the LD30 (upper row) or to the LD80 (lower row). The coloured rectangles in the left panel highlight the different therapy cycles. The dotted black lines in the central and right panels highlight the values of the time averages of the total number of cells and of the mean phenotypic state.

Figure 13 :

 13 Figure 13: Comparing the administration of high drug doses separated by drug-free periods with the continuous administration of a relatively low drug dose. Left panel. Drug delivery schedules under consideration. The blue, green and yellow lines correspond to 12 therapy cycles of 1 week whereby the MTD (i.e. the LD80) is administered either on day 1 [u 1 (t), blue line], day 4 [u 3 (t), green line] or day 7 [u 4 (t), yellow line] of every cycle. The cyan line [i.e. u 2 (t)] refers to the constant infusion for 84 days of a lower amount of drug such that the total dose administered over the whole treatment is the same as that of the other drug delivery schedules. Central panel. Comparison between the average cell population size (in units of 10 8) obtained under the optimal dosing regimen shown in the lower row, left panel of Figure9[i.e. u opt (t)] and the drug delivery schedules u 1 (t), u 2 (t), u 3 (t) and u 4 (t). Right panel. Comparison between the final cell population size (in units of 10 8 ) obtained under the optimal dosing regimen shown in the left panel of Figure10[i.e. u opt (t)] and the drug delivery schedules u 1 (t), u 2 (t), u 3 (t) and u 4 (t).

  ds dz and equation (3.11) follows immediately.

Figure S2 :

 S2 FigureS2: Optimal dosing regimen for 12 therapy cycles of 1 week with 3 days of rest per week with the goal being to minimise the total number of cells at the end of the treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.2), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table1, the values of the terms f 0 and µ 0 given in equation (4.1), ρ 0 = ρ ∞ ū=0 /4 with ρ ∞

  ū=0

Table 1 :

 1

	Parameter	Biological meaning	Value
	γ	Maximum fitness	0.66 day -1
	η	Selection gradient	0.132 day -1
	d	Rate of death due to intrapopulation competition 0.66 × 10 -8 cells -1 day -1
	β	Rate of epimutations	0.001 day -1

Values of the model parameters used to carry out numerical simulations. These values have been estimated based on existing data from in vitro experiments on the phenotypic evolution of HL60 leukemic cells exposed to vincristine

[START_REF] Pisco | Nondarwinian dynamics in therapy-induced cancer drug resistance[END_REF]

.

  If f i (0) = H i then f i (t) = H i for all t ∈ [0, τi ] [cf. the result given by equation (B.8)] and choosing f i (t) = H i in equation (B.12) yieldsµ i (t) = Θ i + (µ i (0) -Θ i ) exp -2 β 1/2 HTherefore, we need to consider further only two cases, namely f i (0) < H i and f i (0) > H i . Before doing this, we rewrite equation (B.12) as Case f i (0) < H i . If f i (0) < H i , using our solution for f i (t) given by equation (B.10) we can write

			β f i (t)		1/2	=	β H i		1/2	coth 2β 1/2 H	1/2 i	t + arctanh (f i (0)/H i ) 1/2
									=	1 2H i	d dt	log sinh 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i ) 1/2	,
	and so for ξ, t ∈ [0, τi ] with ξ < t
		H i	ξ	t	2		β f i (s)	1/2	ds = log	 	sinh 2β 1/2 H i 1/2 sinh 2β 1/2 H 1/2 i ξ + arctanh (f i (0)/H i ) 1/2 t + arctanh (f i (0)/H i ) 1/2	  .
																			1/2
																			u i , for t ∈ (0, τi ].	(B.12)
																			1/2 i	t , for t ∈ [0, τi ].	(B.13)
	d dt	µ i (t) exp H i		0	t	2	β f i (s)	1/2	ds		= 2	β f i (t)	1/2	u i exp H i	0	t	2	β f i (s)	1/2	ds , for t ∈ (0, τi ].
																			(B.14)

  Eliminating the integral inside the exponential in equation (B.14) and multiplying both sides of the resulting equation by K we find that for t ∈ (0, τi ]

								1/2
								ds = log sinh 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i )	1/2	-log (K) .
	d dt	i µ i (t) sinh 2β 1/2 H 1/2	1/2
	= 2	β f i (t)	1/2	u i sinh 2β 1/2 H	1/2 i	t + arctanh (f i (0)/H i ) 1/2
	= 2u i	β H i	1/2	coth 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i ) 1/2	sinh 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i )	1/2
	= 2u i	β H i	1/2	cosh 2β 1/2 H	1/2 i	t + arctanh (f i (0)/H i ) 1/2
	=	u i H i	d dt	sinh 2β 1/2 H	1/2 i	t + arctanh (f i (0)/H i ) 1/2	,
	and so for t ∈ (0, τi ]
						µ i (t) sinh 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i )	1/2	-µ i (0)ζ i
								= Θ i sinh 2β 1/2 H i 1/2	t + arctanh (f i (0)/H i )	1/2	-ζ i ,	(B.15)
	where
									ζ i = sinh arctanh (f i (0)/H i ) 1/2

t + arctanh (f i (0)/H i )

  1/2 H Case f i (0) > H i . If f i (0) > H i ,using our solution for f i (t) given by equation (B.11) we can write

											1/2 i	t + arctanh	f i (0) H i	1/2	, for t ∈ [0, τi ].	(B.16)
		β f i (t)	1/2	=	β H i	1/2	tanh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2
						=	1 2H i	d dt	log cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2	,
	and for ξ, t ∈ [0, τi ] with ξ < t
	H i	ξ	t	2	β f i (s)	1/2	ds = log	 	cosh 2β 1/2 H cosh 2β 1/2 H i 1/2 i 1/2 ξ + arctanh (H i /f i (0)) 1/2 t + arctanh (H i /f i (0)) 1/2	  .

  Eliminating the integral inside the exponential in equation (B.14) and multiplying both sides of the resulting equation by L we find that In the case when f i (0) > H i , substituting the expressions for f i (t) and µ i (t) given, respectively, by equations (3.13) and (3.14) into equation (D.1) yieldsQ i (t) = Γ i -H i λ 2 i (µ i (0) -Θ i ) 2 sech 2 2β 1/2 H-log(λ i ) , for t ∈ [0, τi ].

												1/2 i	t + arctanh (H i /f i (0)) 1/2
										-	1 2	d dt	log cosh 2β 1/2 H	1/2 i	t + arctanh (H i /f i (0)) 1/2	, for t ∈ [0, τi ]
	and so, noting that tanh (z) = sech 2 (z), we obtain
		0	t	Q i (s)ds = Γ i t -	H i 1/2 2β 1/2 (µ i (0) -Θ i ) λ 2 i	2 tanh 2β 1/2 H i 1/2	t + arctanh	H i f i (0)	1/2	-	f i (0) H i	1/2
								-	1 2			log cosh 2β 1/2 H	1/2 i	t + arctanh	f i (0) H i	1/2
												1/2
												ds = log cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2	-log(L).
	d dt			µ i (t) cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2
	= 2		β f i (t)	1/2	u i cosh 2β 1/2 H	1/2 i	t + arctanh (H i /f i (0)) 1/2
	= 2u i	β H i	1/2	tanh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2	cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2
	= 2u i	β H i	1/2	sinh 2β 1/2 H	1/2 i	t + arctanh (H i /f i (0)) 1/2
	=	u i H i	d dt	cosh 2β 1/2 H	1/2 i	t + arctanh (H i /f i (0)) 1/2	,
	and so for t ∈ [0, τi ]
									µ i (t) cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0)) 1/2	-µ i (0)λ i
												=	u i H i	cosh 2β 1/2 H i 1/2	t + arctanh (H i /f i (0))	1/2	-λ i ,	(B.17)
	where
												λ i = cosh arctanh (H i /f i (0)) 1/2	.

* BDH acknowledges support from the Australian Research Council (DP140100339). LA and TL gratefully acknowledge support of the project PICS-CNRS no. 07688 and the French "ANR blanche" project Kibord: ANR-13-BS01-0004.

A Proof of Proposition 1

Substituting the definitions given by equations (2.3) and (2.4) into equation (2.2) and making use of the following definitions

one can rewrite equation (2.2) as

Building upon the results presented by Chisholm et al. [START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to population-level consequences[END_REF] and Lorenzi et al. [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF], we make the ansatz

We note that, since cosh (arctanh(z)) = 1

we have

We can obtain from equation (B.17) an explicit solution for µ i (t) when f i (0) > H i :

Taken together, the results given by equations (B.13), (B. [START_REF] Souza | Chemotherapy dosing schedule influences drug resistance development in ovarian cancer[END_REF]) and (B.18) allow one to obtain the recursive scheme given by equation (3.14).

C Proof of Proposition 3

Under the assumption on u(t) given by equation (3.17 Under the same notation as in the proof of Proposition 2, we have

In the case when f i (0) < H i , substituting the expressions for f i (t) and µ i (t) given, respectively, by equations (3.13) and (3.14) into equation (D.1) yields

and so, noting that coth (z) = -cosech 2 (z), we find

E Supplementary figures

Figure S1: Optimal dosing regimen for 12 therapy cycles of 1 week with the goal being to minimise the total number of cells at the end of the treatment. Dynamics of the total number of cells ρ(t) in units of 10 8 (central panel) and the mean phenotypic state of the cell population µ(t) (right panel) under the optimal drug delivery schedule u(t) (left panel) which is obtained by solving numerically the optimal control problem given by equation (5.6), with J defined by equation (5.2), subject to the dynamic constraints given by the system of differential equations (3.9) and the control constraints given by equation (5.8). Simulations have been carried out over a time window corresponding to 84 days using the parameter values listed in Table 1, the values of the terms f 0 and µ 0 given in equation (4.1), ρ 0 = ρ ∞ ū=0 /4 with ρ ∞ ū=0 given by equation (3.2), and the additional assumptions given by equation (5.10). Focussing on a treatment protocol which consists of 12 therapy cycles of 1 week with a constraint on the total amount of drug that can be administered per cycle, which is the same for all the therapy cycles, we have specified the constraint given by equation (5.8) as in equation (5.11) with the value of C corresponding to the LD80. The coloured rectangles in the left panel highlight the different therapy cycles.