www.biomathforum.org/biomath/index.php/biomath REVIEW ARTICLE A survey of adaptive cell population dynamics models of emergence of drug resistance in cancer, and open questions about evolution and cancer Jean Clairambault∗, Camille Pouchol† ∗INRIA Paris & Sorbonne Université, UMR 7598, LJLL, BC 187, 75252 Paris Cedex 05, France jean.clairambault@ljll.math.upmc.fr †Department of Mathematics KTH, Brinellvägen 8, 114 28 Stockholm, Sweden pouchol@kth.se Received: 9 March 2019, accepted: 14 May 2019, published: 24 May 2019 Abstract—This article is a proceeding sur- vey (deepening a talk given by the first author at the BioMath 2019 International Conference on Mathematical Models and Methods, held in Będlewo, Poland) of mathematical models of cancer and healthy cell population adaptive dynamics exposed to anticancer drugs, to describe how cancer cell populations evolve toward drug resistance. Such mathematical models consist of par- tial differential equations (PDEs) structured in continuous phenotypes coding for the ex- pression of drug resistance genes; they in- volve different functions representing targets for different drugs, cytotoxic and cytostatic, with complementary effects in limiting tu- mour growth. These phenotypes evolve con- tinuously under drug exposure, and their fate governs the evolution of the cell population under treatment. Methods of optimal control are used, taking inevitable emergence of drug resistance into account, to achieve the best strategies to contain the expansion of a tu- mour. This evolutionary point of view, which re- lies on biological observations and resulting modelling assumptions, naturally extends to questioning the very nature of cancer as evo- lutionary disease, seen not only at the short time scale of a human life, but also at the billion year-long time scale of Darwinian evo- lution, from unicellular organisms to evolved multicellular organs such as animals and man. Such questioning, not so recent, but recently revived, in cancer studies, may have conse- quences for understanding and treating can- cer. Copyright: c©2019 Clairambault et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Jean Clairambault, Camille Pouchol, A survey of adaptive cell population dynamics models of emergence of drug resistance in cancer, and open questions about evolution and cancer, Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 1 of 23 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... Some open and challenging questions may thus be (non exhaustively) listed as: • May cancer be defined as a spatially lo- calised loss of coherence between tissues in the same multicellular organism, ‘spa- tially localised’ meaning initially starting from a given organ in the body, but also possibly due to flaws in an individual’s epigenetic landscape such as imperfect control of differentiation genes? • If one assumes that “The genes of cellular cooperation that evolved with multicel- lularity about a billion years ago are the same genes that malfunction in can- cer.” (Davies and Lineweaver, 2011), how can these genes be systematically investi- gated, looking for zones of fragility - that depend on individuals - in the ‘tinker- ing’ (F. Jacob, 1977) evolution is made of, tracking local defaults of coherence? • What is such coherence made of and to what extent is the immune system responsible for it (the self and differen- tiation within the self)? Related to this question of self, what parallelism can be established between the development of multicellularity in different species pro- ceeding from the same origin and the development of the immune system in these different species? Keywords-Cell population dynamics; struc- tured models; Darwinian evolution; drug- induced drug resistance; cancer therapeutics; optimal control I. Introduction: Motivation from and focus on drug resistance in cancer Slow genetic mechanisms of ‘the great evolu- tion’ that has designed multicellular organisms, together with fast reverse evolution on smaller time windows, at the scale of a human disease, may explain transient or established drug re- sistance. This will be developed around the so- called atavistic hypothesis of cancer. Plasticity in cancer cells, i.e., epigenetic [30], [64] (much faster than genetic mutations, and reversible) propension to reversal to a stem- like, de-differentiated phenotypic status, result- ing in fast adaptability of cancer cell popula- tions, makes them amenable to resist abrupt drug insult (high doses of cytotoxic drugs, ion- ising radiations, very low oxygen concentrations in the cellular medium) as response to cellular stress. Intra-tumour heterogeneity with respect to drug resistance potential, meant here to model between-cell phenotypic variability within can- cer cell populations, is a good setting to rep- resent continuous evolution towards drug re- sistance in tumours. This is precisely what is captured by mathematical (PDE) models struc- turing cell populations in relevant phenotypes, relevant here meaning adapted to describe an environmental situation that is susceptible to abrupt changes, such as introduction of a deadly molecule in the environment. Beyond classical (in ecology) viability and fecundity, reversible plasticity for cancer cell populations may also be set as one of such phenotypes. Such structured PDE models have the advan- tage of being compatible with optimal control methods for the theoretical design of optimised therapeutic protocols involving combinations of cytotoxic and cytostatic (and later possibly epi- genetic [89]) treatments. The objective function of such optimisation procedure being chosen as minimising a cancer cell population number, the constraints will consist of minimising unwanted toxicity to healthy cell populations. The innova- tion in this point of view is that success or failure of therapeutic strategies may be evaluated by a mathematical looking glass into the hidden core of the cancer cell population, in its potential of adaptation to cellular stress. The poor understanding of the determinants of drug resistance in cancer at the epigenetic level thus far, and the unexplained failure - or partial failure - of initially promising treatments such as targeted therapies and immunotherapies make it mandatory, from our point of view, to examine cancer, its evolution and its treatment at the level of a whole multicellular organism that - locally, to begin with - progressively lacks Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 2 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... its within- and between-tissue cohesion. II. Biological background A. The many facets of drug resistance in cancer Drug resistance is a phenomenon common to various therapeutic situations in which an external pathogenic agent is proliferating at the expense of the resources of an organism: an- tibiotherapy, virology, parasitology, target pop- ulations are able to develop drug resistance mechanisms (e.g., expression of β-lactamase in bacteria exposed to amoxicillin). In cancer, there is in general no external pathogenic agent (even though one may have favoured the disease) and the target cell populations share much of their genome with the host healthy cell population, making overexpression of natural defence phe- nomena easy (e.g., ABC transporters in cancer cells). Note that drug resistance (and resistance to radiotherapy) is one of the many forms of fast resistance to cellular stress, possibly coded in ‘cold’, i.e., strongly preserved throughout evo- lution, rather than in ‘hot’, i.e., mutation-prone, genes [107]. At the molecular level in a single cell (that is per se insufficient to explain the emergence of drug resistance), overexpression of ABC trans- porters, of drug processing enzymes, decrease of drug cellular influx, etc. [38] are relevant to describe endpoint molecular resistance mecha- nisms. At the cell population level, representing drug resistance by an abstract continuous vari- able x standing for the level of expression of a resistance phenotype (in evolutionary game the- ory [4]: a strategy of the population) is adapted to describe continuous evolution from total sen- sitivity (x = 0) towards total resistance (x = 1). Is such evolution towards drug resistance due to sheer Darwinian selection of the fittest by mutations in differentiation at cell division or, at least partially, due to phenotype adaptation in individual cells? This is by no means clear from biological experiments. In particular, it has been shown in [88] that emergence of drug re- sistance may be totally reversible, and, further- more, that it may be completely dependent on the expression and activity of epigenetic control drugs (DNA methyltransferases). This has been completed by molecular studies of the role of repeated sequences in drug tolerance in [42]. B. Ecology, evolution and cancer in cell popula- tions “Nothing in biology makes senses except in the light of evolution.” (Theodosius Dobzhansky, 1964 [27]) The animal genome (of the host to cancer) is rich and amenable to adaptation scenarios that, especially under deadly environmental stress, may recapitulate salvaging developmental sce- narios - in particular blockade of differentiation or dedifferentiation, allowing better adaptabil- ity but resulting in insufficient cohesion of the ensemble - that have been abandoned in the process of the great evolution [45], [46] from unicellular organisms (aka Protozoa) to coher- ent Metazoa [23] (aka multicellular organisms). In cancer populations, enhanced heterogeneity with enhanced proliferation and poor differenti- ation results in a high phenotypic or genetic di- versity of immature proliferating clonal subpop- ulations, so that drug therapy may be followed, after initial success, by relapse due to selection of one or more resistant clones [25]. As regards ecology and evolution, genetics and epigenetics: ecology is concerned with thriv- ing or dying of living organisms in populations in the context of their trophic environment. Evolution is concerned by the somatic changes, either inscribed by genetic mutations of base pairs in the marble of their DNA, or only - and sometimes geneticists in that case dismiss the term evolution, preferring adaptation - re- versibly (however transmissible to the next gen- erations) modified by silencing or re-expressing genes by means of grafting methyl or acetyl radicals on on the base pairs of the DNA or on the aminoacids that constitute the chromatin (i.e., histones) around which the DNA is wound. In the latter case, such evolution is determined by epigenetic mechanisms, i.e., mechanisms that Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 3 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... do not change the sequence of the base pairs, but only locally modify their transcriptional function. At the level of a same multicellular or- ganism, such so-called epimutations govern the succession of events that physiologically result in cell differentiations. At this stage, one must clearly distinguish three different meanings of the words evolution and differentiation/maturation: 1) at the short- term time scale of a cancer cell population, evolution by (plastic) adaptation of phenotype or by mutations - both phenomena may be encountered - to a changing environment such as infusion of an anti-cancer drug [32]; 2) at the mid-term time scale of a developing animal or human organism, programmed succession of cell differentiations, leading from one original cell to the circa 200-400 different cell types (in the human case) that constitute us as multicellular organisms (development): no mutations, only differentiations (aka maturations), i.e., epige- netic changes within the same genome until each cell reaches its completely physiologically mature state (an example of an ODE model of differentiation may be found, e.g., in [44]); 3) at the very long-term (billions of years) time scale of Darwinian evolution of species, succession of mutations from protozoa until evolved metazoa. Epigenetic changes in a cell population are not rare events and may be fast, operating under environmental pressure by means of epigenetic control enzymes (methylases and acetylases, DNA methyltransferases, etc.), and they are reversible, however likely less quickly than they have occurred [88]. It is also likely, and indeed this has been shown in some cases, that follow- ing such reversible epigenetic events, rare events that are mutations (not so rare in the context of genetic instability that often characterises cancer cell divisions) stochastically happen, fix- ing in the DNA an acquired advantage in the context of a changing environment. Conversely, mutations in parts of the genome that code for epigenetic control enzymes may determine epi- genetic changes, metaphorically representable in the Waddington epigenetic landscape [45], [46], [106]. Such genetic to epigenetic modifications and vice versa are discussed in [14], [36], [37]. Also note that the relationships between ecol- ogy, evolution and cancer are extensively devel- oped in the book [102]. C. The atavistic theory of cancer There has been some debate in the past 20 years about two opposed views of cancer, the classic one, advocated by P. Nowell [71] states that cancer starts from a single “renegade” cell (a cheater), that by a succession of muta- tions initiates cancer (Somatic Mutation The- ory, SMT), being followed by strict Darwinian selection of the fittest, while the less admitted Tissue Organisational Field Theory (TOFT), advocated in particular by A. Soto and C. Son- nenschein [90], contends that cancer is a matter of deregulated ecosystem, amenable to eradica- tion by changing the tumour ecosystem. The atavistic theory of cancer is completely different from those two, in as much as it relates cancer to a regression in the billion year-long evolu- tion of multicellularity. The idea that cancer is a form of backward evolution from organised multicellularity toward unicellularity, stalled at the poorly organised forms of multicellularity tumours consist of (as cancer cell populations, escaping the collective control present in delicate organismic organisations that constitute coher- ent multicellularity, continuously reinvents the wheel of multicellularity, starting from scratch for their own sake) is not new and it has at least been proposed by T. Boveri in 1929 [7] and L. Israel in 1996 [47]. However, it has regained visibility thanks to the documented and simultaneous studies by physicists P.C.W. Davies and C.H. Lineweaver [23], [55] and oncol- ogist M. Vincent [103], [104], followed by various subsequent studies, constituting a new body of knowledge [11], [19], [97], [98], [108] under the name of atavistic theory of cancer. This theory, or hypothesis, postulates that, although cancer reinvents the wheel of multicellularity, it has at its disposal for this task “an ancient Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 4 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... toolkit of pre-existing adaptations” that makes it fundamentally differ from classical Darwinian evolution [23]. In this respect, cancer is clearly “more an archeoplasm than a neoplasm” (Mark Vincent [103]). What is relatively new in this theory, com- pared with the previously cited ones, is the idea that an intermediate set of coarse forms of multicellularity (which they call “Metazoa 1.0”), lacking coherent control of intercellular - and between cell populations - cooperativ- ity and proliferation, qualified a “robust toolkit for the survival , maintenance and propagation of non-differentiated or weakly-differented cells” is a safety state to which “Metazoa 2.0” (us, in particular) revert when our sophisticated form of multicellularity goes astray in cancer. Such events are due to failures in control of evolved cooperativity genes, and this incoher- ent, chaotic, poorly organised “Metazoa 1.0” system endows the tissues in which it is installed with high phenotype adaptability, aka cancer plasticity, on which tumour development relies. Such plasticity makes tumour cells in particular able to exploit for their own sake, plasticity resulting in resistance to cytotoxic drugs, epige- netic enzymes that were physiologically designed to control finely tuned cell differentiations, in acquired resistance to cancer treatments. This illuminating view of cancer, according to which the genes that malfunction are precisely “the genes of cellular cooperation that evolved with multicellularity about a billion years ago” (Paul Davies and Charles Lineweaver [23]), and the reason of resistance is to be found in ancient, well preserved, genes of our DNA, is however quite often not admitted by many biologists of cancer who strongly believe in the strictly Darwinian nature of evolution in cancer cell populations [33], [34], [39], [40], without any kind of such “genomic memory”. Compatible with the atavistic hypothesis that postulates such backward evolution, a possible scenario suggests that cancer may start with a local deconstruction of the epigenetic con- trol of cell differentiation (that is an essential piece of the coherence of multicellularity, e.g., in haematopoiesis, by genes TET2, DNMT3A, ASXL1), followed by deregulation of cooper- ativity between cell populations (essential to division of work in a multicellular organism) initiated by disruption of transcription factors responsible for differentiation (e.g., by genes RUNX1, CEBPα, NPM1) and finally deregula- tion of the determinants of the strongest and most ancient bases of multicellularity, prolif- eration and apoptosis (e.g., by genes FLT3, KIT, and genes of the RAS pathway). Even though many cancer biologists are reluctant to endorse this scenario, biological observations ex- ist, showing that a scenario of successions of mutations may be found in fresh blood sam- ples of patients with acute myeloid leukaemia, phylogenetically recapitulating such hierarchi- cally ordered deconstruction of the multicellu- lar haematopoietic structure, from the finest (epigenetic) to the coarsest (proliferation and apoptosis) elements of the construction of mul- ticellularity [43]. From the point of view of therapeutic appli- cations, the atavistic theory of cancer has the consequence that, even though those genes of cooperativity that are altered in cancer (the “multicellularity gene toolkit of Metazoa 1.0”) have taken one billion years of Darwinian evo- lution to achieve (by ‘tinkering’ [49]) coher- ent evolved multicellular organisms, they are nevertheless in finite number and can be sys- tematically investigated, as has been initiated in phylostratigraphic studies led by Tomislav Domazet-Lošo and Diethard Tautz [28], [29]. Such systematic between-species phylogenetic biocomputer studies should open observation windows onto altered genes in patients and their possible correction in the future. From the point of view of mathematical mod- elling, the fact that ancient genes of survival have been developed in the course of evolution to make individual cells, and later coherently heterogeneous and nevertheless communicating Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 5 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... together (failures in intercellular communica- tions [99], [100], [101] incidentally being also a possible source of default of cooperativity in cell populations), and may be conserved as silent capacities in our genome, only waiting to be unmasked by epigenetic enzymes [88] put on the service of survival in highly plastic cancer cells [107], gives reasonable biological support to the notion of cell populations structured in phenotype of survival and of drug resistance. How such ancient genes (‘cold genes’ [107]) have been preserved in our genome while serving in rare and extreme environmental conditions only is not clear (in principle, genes that are not expressed are prone to disappear), however observations reported in [107] propose that an- cient genes evolve more slowly than younger ones. Hence, preserved in the genomic memory as survival genes, revivable in plastic cancer cell populations (plastic here meaning that they have easy access to epigenetic enzymes to change their phenotype under environmental pressure), their level of expression may offer a basis for evolvability and reversibility, under environmen- tal pressure, of continuous phenotypes structur- ing the heterogeneity (aka biological variability) of cancer cell populations that is developed in mathematical models of adaptive cell population dynamics. III. Models of adaptive dynamics A. Models structured in resistance phenotype The simplest model of a resistance phenotype- structured cell population may be described by a nonlocal Lotka-Volterra-like integro-differential equation, here x ∈ [0, 1] representing a contin- uous resistance phenotype, from x = 0, total sensitivity to the drug, until x = 1, total insen- sitivity: ∂n ∂t (t,x) = ( r(x) −d(x)ρ(t) ) n(t,x), with ρ(t) := ∫ 1 0 n(t,x) dx and n(0,x) = n0(x). Note that this simple integro-differential equation may, when this makes biological sense, be gen- eralised to a reaction-diffusion-advection (RDA) one written as ∂n ∂t (t,x) + ∂ ∂x {v(x)n(t,x)} = β ∂2 ∂x2 n(t,x) + {r(x) −d(x)ρ(t)}n(t,x). We assume reasonable (L∞) hypotheses on r and d, and n0 ∈ L1([0, 1]). Phenotype- dependent functions r and d stand for intrinsic proliferation rate and intrinsic death rate due to within-population competition for space and nutrients, respectively. Note that space is repre- sented here only in the abstract nonlocal logistic term d(x)ρ(t). It is nevertheless possible to mix phenotype and actual Cartesian space variables to structure the population, as will be shown later. One can then prove for the simple integro- differential model the asymptotic behaviour the- orem: Theorem 1. [24], [48], [74] (i) ρ converges to ρ∞, the smallest value ρ such that r(x) − d(x)ρ ≤ 0 on [0, 1] (i.e., ρ∞ = max[0,1] rd). (ii) The population n(t, ·) concentrates on the phenotype set { x ∈ [0, 1], r(x) −d(x)ρ∞ = 0 } . (iii) Furthermore, if this set is reduced to a singleton x∞, then n(t, ·) ⇀ ρ∞δx∞ in M1(0, 1). (the measure space M1(0, 1) being the dual for the supremum norm of the space of continu- ous real-valued fonctions on [0, 1]; note the Dirac mass on the RHS, convergence is here meant in the sense of measures.) Although in the one-population case, as stated above, a direct proof of convergence based on proving that ρ(t) is BV on the half-line, from which concentration easily follows from exponential growth, it is interesting to note, as this argument can be used in the case of two interacting populations, that a global proof based on the design of a Lyapunov function gives at the same time convergence and concentration: Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 6 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... choosing any measure n∞ on [0, 1] with support in argmax r d such that∫ 1 0 n∞(x) dx = ρ∞ = max [0,1] r d , and for an appropriate weight w(x) (in fact 1 d(x) )setting V(t)= ∫ 1 0 w(x){n(t,x)−n∞(x)−n∞(x) ln n(t,x)}dx, one can show that dV dt = −(ρ(t) −ρ∞)2 + ∫ 1 0 w(x){r(x)−d(x)ρ∞}n(t,x) dx, which is always nonpositive, tends to zero for t → ∞, thus making V a Lyapunov function, and showing at the same time convergence and concentration. Indeed, in this expression, the two terms are nonpositive and their sum tends to zero; the zero limit of the first one accounts for convergence of ρ(t), and the zero limit of the second one accounts for concentration in x (on a zero- measure set) of lim t→+∞ n(t,x). Starting from this simple model, one can gen- eralise it to the case of two interacting cell pop- ulations, cancer (nC) and healthy (nH), again using a nonlocal Lotka-Volterra setting, with two different drugs, u1, cytotoxic (= cell-killing drug, towards which resistance evolves according to the continuous phenotype x ∈ [0, 1]) and u1, cytostatic (only thwarting proliferation without killing cells): ∂ ∂t nH(t,x) =[ rH(x) 1+kHu2(t) −dH(x)IH(t)−u1(t)µH(x) ] nH(t,x), ∂ ∂t nC(t,x) =[ rC(x) 1+kCu2(t) −dC(x)IC(t)−u1(t)µC(x) ] nC(t,x). (1) The environment in the logistic terms is defined by: IH(t) = aHH.ρH(t) + aHC.ρC(t), IC(t) = aCH.ρH(t) + aCC.ρC(t), with aHH > aHC, aCC > aCH (higher within- species than between-species competition) and ρH(t) = ∫ 1 0 nH(t,x) dx, ρC(t) = ∫ 1 0 nC(t,x) dx. The cytotoxic drug terms, tuned by drug sen- sitivity functions µC and µH, act as added death terms to the logistic term, whereas the cyto- static drug terms act by inhibiting the intrinsic proliferation rates rC and rH. Functions µC and µH obviously have to be decreasing functions of x, and so, less obviously, but representing a trade-off between survival and proliferation (“cost of resistance”), have to be rC and rH. As regards dC and dH, no modelling choice imposes itself; however, in order to make the function r d globally decreasing and thus, in the absence of drug, obtain its maximum around zero, it was assumed in this study that it is a non- decreasing function of x. Biologically, this means that the more resistant a cell is, the stronger opposition to its proliferation it encounters in its own species, cancer or healthy, which is another way, coherent with the modelling choice made on rC and rH, to express a cost of resistance. In this 2-population case, following an ar- gument by Pierre-Emmanuel Jabin and Gaël Raoul [48], one can prove, as in the 1-population case, at the same time convergence and concen- tration by using a Lyapunov functional of the form∫ w(x){n(t,x) −n∞(x) −n∞(x) ln n(t,x)} dx. We have also in this case the asymptotic behaviour theorem: Theorem 2. [81], [83] Assume that u1 and u2 are constant: u1 ≡ ū1, and u2 ≡ ū2. Then, for any positive initial population of healthy and of tumour cells, (ρH(t),ρC(t)) converges to the Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 7 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... equilibrium point (ρ∞H ,ρ∞C ), which can be exactly computed as follows: Let a1 ≥ 0 and a2 ≥ 0 be the smallest nonnega- tive real numbers such that rH(x) 1 + αHū2 − ū1µH(x) ≤ dH(x)a1 and rC(x) 1 + αCū2 − ū1µC(x) ≤ dC(x)a2. Then (ρ∞H ,ρ∞C ) is the unique solution of the invertible (aHH.aCC > aCH.aHC) system I∞H = aHHρ ∞ H + aHCρ ∞ C = a1, I∞C = aCHρ ∞ H + aCCρ ∞ C = a2. Let AH ⊂ [0, 1] (resp., AC ⊂ [0, 1]) be the set of all points x ∈ [0, 1] such that equality holds in the inequalities above. Then the supports of the probability measures νH(t) = nH(t,x) ρH(t) dx and νC(t) = nC(t,x) ρC(t) dx converge respectively to AH and AC as t tends to +∞. In [81], this result is complemented with nu- merical simulations which show the failure of constant administration of high doses of both drugs. The theorem explains the phenomenon: such a strategy makes the cancer cell density concentration on a very resistant phenotype near x = 1. Once most of the mass is close to x = 1, further treatment is hopeless as the tu- mour has become mostly resistant, and it starts increasing again after having first decreased. This can be interpreted as relapse. Note that numerical studies based on a sim- ilar model of adaptive dynamics in a reaction- diffusion version, dealing with the question of relapse, can be found in [16], [17], [56], [57], see also [74], [75] for more theoretical consid- erations. This result extends to two competitively in- teracting populations the result of convergence and concentration for nonlocal Lotka-Volterra phenotype-structured models previously pub- lished in [24], [48], [74]. Note that it assumes the invertibility of the square matrix [aij] where i,j ∈{H,C}. A natural question then arises: is it possible to extend this result to N > 2 interacting populations? This is the object of the study [82], in which the following N-dimensional nonlocal Lotka-Volterra system is set, for which one can look for coexistence of positive steady states (i.e., persistence of all species): ∂ ∂t ni(t,x) =  ri(x) + di(x) N∑ j=1 aijρj(t)  ni(t,x), in which x stands for all xi ∈ Xi for simplicity, each Xi being a compact subset of some Rpi, ri,di smooth enough, and as usual ρi(t) = ∫ Xi ni(t,x) dx. This system generalises to a nonlocal setting classical Lotka-Volterra models (for 2 popula- tions in an ODE setting, see, e.g., Britton [8] or Murray [69]) with ecological cases: mutualistic if aij > 0 and aji > 0, competitive if aij < 0 and aji < 0 , predator-prey-like if aijaji < 0, for the interaction matrix A = [aij] In [82], to which the reader is sent for more details, it is proved that a coexistent positive steady state ρ∞ = [ρ∞1 , . . . ,ρ ∞ i , . . . ,ρ ∞ N ] t exists in RN if and only if, setting I∞ = [I∞1 , . . . ,I ∞ i , . . . ,I ∞ N ] t, where each I∞i = max x∈Xi ri(x) di(x) , the equation Aρ + I∞ = 0 has a solution ρ∞ ∈ RN. Then, under some precise conditions on A, it can be proved, again using the same kind of Lyapunov function as in [81], that the solution to the N-dimensional nonlocal Lotka-Volterra system exists and is globally defined; furthermore, the solution ρ∞ to the equation Aρ + I∞ = 0 in RN is then unique and globally asymptotically stable. As in the 1- and 2-dimensional cases, a result of concentration in phenotype follows, with moreover an estimation of the speeds of convergence and of concentration. Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 8 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... B. Modelling mutualistic tumour-stroma inter- actions Noting that breast and prostate tumours are accompanied in their stroma by so-called cancer-associated adipocytes (CAAs) or cancer- associated fibroblasts (CAFs) [18], [26], [60], which favour cancer growth, likely by exchang- ing bidirectional messenger molecules, one can model such mutualistic interactions in a nonlo- cal Lotka-Volterra way: ∂ ∂t nA(t,x) =[rA(x)−dAρA(t)+sA(x)ϕC(t)]nA(t,x) ∂ ∂t nC(t,y) =[rC(y)−dCρC(t)+sC(y)ϕA(t)]nC(t,y) with ρA(t) = ∫ nA(t,x) dx, ρC(t) = ∫ nC(t,y) dy, ϕA(t) = ∫ ψA(x)nA(t,x) dx and ϕC(t) = ∫ ψC(y)nC(t,y) dy, for some weight functions ψA and ψC (that in the absence of known data may be chosen as simply affine), and some given initial conditions nA(0,x) = n0A(x),nC(0,y) = n0C(y) for all (x,y) in [0, 1]2, x standing for transformation towards a CAA or CAF state in the adipocyte popula- tion and y standing for strength of malignancy in the cancer cell population. This model is studied, theoretically and numerically, in [80] in its generalised reaction-diffusion-advection form (see above) with explicit functions and initial functions n0A,n0C assumed to be Gaussian. In a setting in which mutualistic interactions between two cell species, one of them being ini- tially healthy, but susceptible to become cancer- ous, namely proliferating haematopoietic stem cells and early progenitors nh, in the mandatory presence of the other species ns, supporting stromal cells, a model closely related to the previous one is presented [70]. It writes  ∂tnh(t,x) = [ rh(x)−ρh(t)−ρs(t)+α(x)Σs(t) ] nh, ∂tns(t,y) = [ rs(y)−ρh(t)−ρs(t)+β(y)Σh(t) ] ns. This system is completed with initial data nh(0,x) = nh0(x) ≥ 0, ns(0,y) = ns0(y) ≥ 0. Here the assumptions and notations are • ρh(t) := ∫ b a nh(t,x)dx, ρs(t) := ∫ d c ns(t,y)dy are the total populations of HSCs and their supporting stromal cells MSCs, re- spectively, x representing a malignancy po- tential in haematopoietic cells and y a trophic potential in stromal cells. • The functions Σh(t) := ∫ b a ψh(x)nh(t,x) dx, and Σs(t) := ∫ d c ψs(y)ns(t,y) dy denote an assumed chemical signal (Σh) from the hematopoietic immature stem cells (haematopoietic stem cells, HSCs) to their supporting stroma (mesenchymal stem cells, MSCs), i.e., “call for support” and conversely, a trophic message (Σs) from MSCs to HSCs. The weight functions ψh,ψs are nonnegative and defined on (a,b) and (c,d), intervals of the real line. • The function rh ≥ 0 represents the intrin- sic (i.e., without contribution from trophic messages from MSCs) proliferation rate of HSCs. Assume that rh is non-decreasing, rh(a) = 0 and rh(b) > 0. • The function α ≥ 0, satisfying α′ ≤ 0 and α(b) = 0, is the sensitivity of HSCs to the trophic messages from supporting cells. • For the function rs ≥ 0, it is assumed that r′s(y) ≤ 0. The function β(y) ≥ 0 with β′(y) ≥ 0 represents the sensitivity of the stromal cells MSCs to the (call for support) message coming from HSCs. Some examples for rh, α are given by rh = r∗h(x−a) or rh = r ∗ h(x−a) 2, α(x) = α∗(b−x) with positive constants r∗h,α ∗, ψs(y) = y and ψh(x) = x. The reader is sent to [70] for a detailed study of this model. In particular, theoretical condi- tions for extinction, invasion or possible stable Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 9 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... coexistence of a leukaemic clone emerging in an initial healthy HSC population together with a maintained healthy fraction of it, with numerical simulations, are given in this study. They are related to convexity or concavity properties of functions of the model describing proliferation of the population, rh and α, the same kind of evolution being possible in the stromal cell population. C. Models structured in phenotype and space Although purely space-structured models lack the necessary heterogeneity in phenotype to take into account continuous evolution towards drug-induced drug resistance, purely phenotype- structured models lack the possibility to exam- ine possible heterogeneities due to extension of tumours in Cartesian space, in particular due to diffusion of molecules (anticancer drugs and nutrients) in the medium. Hence, provided that something is known of the geometry of the space occupied by a cancer cell population, and this is indeed the case with initial tumours that spontaneously thrive in spheroids, mixing space with phenotype to structure a model of a cancer cell population under drug exposure, to study its behaviour with respect to drug resistance, is a natural way to proceed. Relying on modelling principles developed in [52], [58], integrated in spheroid-like space, such a model is studied in [59]: ∂tn(t,r,x) = [ p(x) 1+µ2c2(t,r) s(t,r) −d(x)%(t,r)−µ1(x)c1(t,r) ] n(t,r,x), −σs∆s(t,r)+ [ γs+ ∫ 1 0 p(x)n(t,r,x)dx ] s(t,r) = 0, −σc∆c1(t,r)+ [ γc+ ∫ 1 0 µ1(x)n(t,r,x)dx ] c1(t,r) = 0, −σc∆c2(t,r)+ [ γc+µ2 ∫ 1 0 n(t,r,x)dx ] c2(t,r) = 0, with zero Neumann conditions at r = 0 (spheroid centre) coming from radial symmetry and Dirichlet boundary conditions at r = 1 (spheroid rim): s(t,r = 1) = s1, ∂rs(t,r = 0) = 0, c1,2(t,r = 1) = C1,2(t), ∂rc1,2(t,r = 0) = 0, where: • The function p(x) is the intrinsic (i.e., in- dependently of cell death) proliferation rate of cells expressing resistance level x due to the consumption of resources. The factor 1 1 + µ2c2(t,r) mimics the effects of cytostatic drugs, which act by slowing down cellular proliferation, rather than by killing cells. The parameter µ2 models the average cell sensitivity to these drugs. • The function d(x) models the death rate of cells with resistance level x due to the competi- tion for space and resources with the other cells. • The function µ1(x) denotes the destruction rate of cells due to the consumption of cytotoxic drugs, whose effects are here summed up directly on mortality. • Parameters σs and σc model, respectively, the diffusion constants of nutrients and cyto- toxic/cytostatic drugs. • Parameters γs and γc represent the decay rate of nutrients and cytotoxic/cytostatic drugs, respectively. The model can be recast in the equivalent form ∂tn(t,r,x) =R ( x,%(t,r),c1,2(t,r),s(t,r) ) n(t,r,x), in order to highlight the role played by the net growth rate of cancer cells, which is described by R ( x,%(t,r),c1,2(t,r),s(t,r) ) := p(x) 1 + µ2c2(t,r) s(t,r) −d(x)%(t,r) −µ1(x)c1(t,r). The following considerations and hypotheses are assumed to hold: • With the aim of translating into mathematical terms the idea that expressing cytotoxic re- sistant phenotype implies resource reallocation Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 10 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... (‘cost of resistance’, i.e., redistribution of ener- getic resources from proliferation-oriented tasks toward development and maintenance of drug resistance mechanism, such as higher expression or activity of ABC transporters in individual cells), p is assumed to be decreasing p(·) > 0, p′(·) < 0. As regards function d, one can note that in this study [59], the advocated modelling choice (d′(·) < 0) is the opposite of the one that was made in [58], a study nevertheless published by the same authors. In [58], the underlying biological reason is possibly that ‘the more re- sistant a cell is, the stronger opposition to its proliferation it encounters in its own species, cancer or healthy, which is another way, coherent with the modelling choice made on rC and rH, to express a cost of resistance’ (see this argu- ment developed above in Subsection III-A). As a matter of fact, the simulations shown in this study [59] always use a constant value for d, and, contrary to [58], no theorem is proposed to the reader of [59], which should induce to actually choose d′(·) ≥ 0. • The effects of resistance to cytotoxic therapies are modeled by the obvious condition that the drug sensitivity function µ1 is non-increasing: µ1(·) > 0, µ′1(·) ≤ 0. The interesting results of this model consist of simulations, illustrated by figures to which the interested reader is referred. D. Models structured in cell-functional variables A puzzling observation on an in-vitro ag- gressive cancer cell culture (PC9, a variant of NSCLC cells) exposed to high doses of anti- cancer drug, experiment reported in [88], is that: 1) even though 99.7% of cells quickly die when exposed to the drug, sparse and tiny subpopula- tions (0.3%) survive, named drug-tolerant per- sisters (DTPs), and for some time just survive, exposed to the same very high concentration of drug; 2) after some time (not precisely defined in the paper), these surviving cells change their phenotype as expressed by membrane markers, and proliferate again, then named drug-tolerant expanded persisters (DTEPs), unabashed in the maintained high drug dose; 3) when the drug is washed out from the cell culture, the cell population reverts to initial drug sensitivity, and such resensitisation occurs ten times more slowly at the DTEP stage than at the DTP stage; 4) if the cell culture is exposed to an inhibitor of the epigenetic enzyme KDM5A together with the drug, be it at the DTP or DTEP stage, DTPs - or DTEPs - die. Such clearly epigenetic and completely reversible mode of resistance, developed in two stages, called for designing a cell population dynamic model structured, not as previously, monotonically in drug resistance gene expression level, but in phenotypes linked to the cell fate, which in cell populations al- ways may be reduced to proliferation, death or differentiation (senescence being a version of delayed death). In the modelling and numerical study [13], 2 phenotypes are thus chosen to take into account the cell population heterogeneity relevant for the experiment: survival potential under extreme environmental conditions (called by ecology theoreticians viability), x, and pro- liferation potential (called fecundity), y. The resulting model is described by the reaction- diffusion-advection equation, that describes the behaviour of a very plastic cell population under exposure to a high dose of anticancer drug: ∂n ∂t (x,y,t) + ∂ ∂y ( v(x,c(t); v̄)n(x,y,t) ) ︸ ︷︷ ︸ Stress-induced adaptation of the proliferation level = β∆n(x,y,t)︸ ︷︷ ︸ Non-genetic phenotype instability + [ p(x,y,%(t))−d(x,c(t)) ] n(x,y,t)︸ ︷︷ ︸ Non local Lotka-Volterra selection • %(t)= ∫ 1 0 ∫ 1 0 n(x,y,t) dx dy, p(x,y,%(t))=(a1+a2y+a3(1−x)) ( 1− %(t) K ) and d(x,c) = c(b1 + b2(1 −x)) + b3 Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 11 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... • The global population term %(t) = ∫ 1 0 ∫ 1 0 n(x,y,t) dx dy occurs in p as a logistic environment lim- iting term (availability of space and nutri- ents). • The drift term w.r.t. proliferation po- tential y represents possible (if v 6= 0) ‘Lamarckian-like’, epigenetic and re- versible, adaptation from PC9s to DTPs; switching from v ≥ 0 to v = 0 here means switching from a possible adaptation scenario to a strictly Darwinian one (it is biologically impossible to decide between the two scenarios). • v(x,c(t); v̄) = −v̄c(t)H(x∗ − x) where t 7→ c(t) is the drug infusion function, and x∗ is a fixed viability threshold. • No-flux boundary conditions. Of note, another, individual-based, model (IBM) yielding the same simulation results (no theorem) is proposed in a complementary way to the interested reader, sent to [13]. The simulation results firstly show total reversibility to drug sensitivity when the drug is withdrawn, and also allow to study the evolution of the two phenotypes in the absence of drug, under drug exposure, and when the drug is withdrawn. Furthermore, the model was put at stake by asking 3 questions: Q1. Is non-genetic instability (Laplacian term) crucial for the emergence of DTEPs? Q2. What can we expect if the drug dose is low? Q3. Could genetic mutations, i.e., an integral term involving a kernel with small support, to replace both adapted drift (advection) and non-genetic instability (diffusion), yield similar dynamics? Consider c(·) = constant and two scenarios: (i) (‘Lamarckian’ scenario (A): the outlaw) Only PC9s initially, adaptation present (v 6= 0) (ii) (‘Darwinian’ scenario (B): the dogma) PC9s and few DTPs initially, no adaptation (v = 0) To make a long story short [13], • Q1. Always yes! Whatever the scenario. • Q2. Low doses result in DTEPs, but no DTPs. • Q3. Never! Whatever the scenario. Can such cell-functional models be used to actually manage drug resistance in the clinic? An idea would be to counter the plastic adap- tation that cancer cell populations show in the presence of high doses of drugs by infusing at the same time as cytotoxic drugs inhibitors of epigenetic enzymes such as KDM5A in [88]. However, even though epigenetic drugs are the object of active research in the pharmaceutic in- dustry [89], the importance of epigenetic control of physiological processes (all differentiation is epigenetic!) and the role of impaired epigenetic controls in impaired cell differentiation, which is a characteristic of cancer, has been stressed [31] makes them delicate to manage in the clinic so far. IV. Optimisation and optimal control A. ODE models and their optimal control in cancer To go beyond the administration of constant doses, one is led to let drug infusion rates vary in time and try to find the best such rates to minimise a given criterion, such as the number of cancer cells at the end of a given time-window. This is the purpose of the mathematical field of optimisation (see, e.g., in another framework [5], [20]) and optimal control, with all its available theoretical and numerical tools. Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 12 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... At this stage, it is noteworthy that the dis- cretisation of the phenotype-structured PDE models introduced so far leads to ODEs, usually of Lotka-Volterra type. To illustrate the idea, let us go back to the prototype integro-differential model ∂n ∂t (t,x) = ( r(x) −d(x)ρ(t) ) n(t,x). If one discretises the phenotype space into Nx + 1 equidistant phenotypes through xi = i∆x, ∆x = 1 Nx the above equation is approxi- mated by the ODE system y ′ i(t) = ( r(xi) −d(xi)ρ(t) ) yi(t), i = 0, . . . ,N where yi(t) ≈ ∆xn(t,xi) and ρ(t) = ∑N j=0 yi(t). This remark is general and applies to the numerical simulation of phenotype-structured PDE models (this is nothing but a semi- discretisation of the corresponding PDE). This point of view also makes the link between ODE models where resistance is represented by a binary variable, or more generally, a dis- crete variable. With a coarser discretisation, the ODE model has few equations and is more amenable to parameter identification, quick nu- merical simulation, but is also less accurate in representing resistance. When it comes to optimal control, ODE mod- els with a discrete representation of resistance have long been studied, either theoretically or numerically [91], [22], [50]. This is one aspect of the rich literature on optimal control for cancer modelling, see the reference book [86]. Note that these ODE models can be made richer, as they may additionally model healthy cells, cells in different compartments of the cell cycle, immune cells, etc. Independently of the number of equations, the investigation of the optimal control problem typically leads to optimal strategies being the concatenation of bang-bang and singular arcs. Bang-bang arcs correspond to drugs being ei- ther given at the maximum tolerated dose or not at all, whereas singular arcs correspond to intermediate doses which can be computed in feedback form from the yi’s. These results are obtained either numerically, or theoretically by applying the Pontryagin Maximum Principle, possibly with higher order criteria (such as the Legendre-Clebsch criterion) and/or with geo- metric optimal control techniques. The usual clinical practice is to use maximum tolerated doses, a strategy which has been called into question as it can lead to an initial drop in tumour size before regrowth due to acquired resistance [25]. This corresponds to bang-bang controls. Instead, alternative and more recent strategies advocate for the infusion of interme- diate doses [35], [73], [93]. Thus, as explained in [53], understanding whether the optimal controls do contain sin- gular arcs is of paramount importance, and might depend from the parameters governing the cost. The recent work [12], where resistance is modelled to be binary, also features parameter regions leading to singular arcs which follow a first arc with maximum tolerated doses. This naturally poses the question of opti- mal scheduling for PDE models of resistance. The corresponding optimal control problems are then significantly harder to solve. Numerically, this is because a fine discretisation leads to computationally demanding algorithms. Theo- retically, as is already the case for a high dimen- sional ODE, it becomes more difficult to obtain precise results on the optimal control strategy, even from an (infinite-dimensional) Pontryagin Maximum Principle. B. Optimal control of phenotype-structured PDE models These difficulties might explain why there are up to date few optimal control results on phenotype-structured PDEs for resistance. Most studies are restricted to constant doses and the optimisation is then performed on the resulting scalar parameters. This can be done by nu- merical investigation of the parameter space as in [16], [58] or theoretically for cases in which explicit solutions are available [3]. Other non- constant infusion strategies mimicking popular Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 13 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... protocols are sometimes also tested in the afore- mentioned works. Finally, in [3], restriction to Gaussian solutions allows the authors to reduce the PDE to a system of three ODEs (for the total mass, the mean and the standard devia- tion). All these studies conclude that the con- tinuous administration of maximum tolerated doses might lead to relapse, and that alternative strategies with lower infusion of drugs might be preferable. There is, up to our knowledge, very little work in the direction of tackling a full optimal control problem for the phenotype-structured PDEs models, without such simplifications as above. The two works we are aware of are [72] and [81], both concerned with the model (1) (see Section III-A). In [72], the model is more complex since genetic instability is introduced in the PDE, modelled by diffusion terms. The goal is to minimise the total number of cancer cells ρC(T), and the overall model is complemented with the constraints • maximum tolerated doses: 0 ≤ u1(t) ≤ umax1 , 0 ≤ u2(t) ≤ u max 2 , • control of the tumour size: ρH(t) ρH(t) + ρC(t) ≥ θHC, (2) • control of the toxic side-effects: ρH(t) ≥ θHρH(0), (3) where 0 < θHC,θH < 1. In order to solve the problem numerically, the approach consists in discretising the whole problem in phenotype and time, thus using a so-called direct method in numerical optimal control [96]. This is equivalent to discretising in time an ODE system which has as many equa- tions as there are discretised phenotypes. The optimal control problem then becomes a high finite-dimensional optimisation problem, which can be handled, for example, by interior point methods. As is common to most numerical optimisation problems, the biggest difficulty lies in choosing the initial guess for the algorithm. The approach of [81] is to solve the optimisation problem with a very coarse discretisation (few unknowns) before scaling the problem up progressively to a fine discretisation. For the generalised model with mutations, such a strategy fails because of the computational cost of Laplacians. To circumvent this, the numerical strategy introduced in [72] is to simplify the PDEs by setting some coefficients to zero, so that the resulting optimal control problem can be solved by a Pontryagin Maximum Principle. Although this problem is non-realistic from the applica- tive a point of view, it provides an excellent starting point for a homotopy procedure which allows to go all the way back to the original more complicated problem, with a very accurate discretisation. An optimal strategy clearly emerges from these two works, when the initial tumour is heterogeneous (as a result of a first standard administration of cytotoxic drugs). The idea is to let the tumour density evolve to a sen- sitive phenotype by using no cytotoxic drugs and intermediate (constant) doses of cytostatic drugs for a long phase, during which the con- straint on tumour size saturates. Only then one takes profit of a sensitive tumour by using the maximum tolerated doses, up until the side- effects constraint saturates. It is then possible to further reduce the tumour size by lowering the cytotoxic dosage. The asymptotic analysis comes in handy in understanding the optimality of such a strategy: the first long phase leads to the convergence of the cancer cell density onto a Dirac mass located on a sensitive phenotype (or a smoothed version of such a Dirac mass when there is a diffusion term). This property allows the authors of [81] to perform a theoretical study of the optimal control problem in a reduced control set where the controls are forced to take constant values during a first long phase. The strategy obtained numerically is then proved to be optimal in a theorem, informally given below. Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 14 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... Theorem 3. [81] When the final time T is large, the optimal solution is such that 1) at the end of the first phase, the density of cancer cells has concentrated on a sensitive phenotype, 2) the optimal strategy is then the concatenation of three arcs • an arc with saturation of the constraint on ρH ρH +ρC . • a free arc with maximum tolerated doses, namely u1 = umax1 and u2 = umax2 , • an arc with saturation of the constraint on ρH and u2 = umax2 . We insist that this result is proved only in the absence of diffusion. The proof relies on the fact that Dirac mass concentration at the end of the first phase allows to replace the PDE system by an 2x2 ODE system, up to an error becoming arbitrarily small as the length of the phase increases. The resulting optimal control problem can then be handled with a Pontryagin Maximum Principle with state constraints. C. Future prospects in optimal control Applying the strategy advocated in [72], [81] requires thinking it in a quasi-periodic manner, and as a strategy relevant after the traditional admnistration of the first dose, which usually induces resistance. The idea would then be to alternate between: 1. a long phase with cytostatic doses and no cytotoxic doses (a drug holiday) to resensi- tise the tumour, 2. a short phase with maximum tolerated doses until the toxicity is considered to have reached its limit, with a possible subsequent switch in dose for the cytotoxic drugs to keep diminishing the tumour size. Such a protocol requires to assess the level of resistance in order to decide when to switch from 1. to 2. and determine when damage to healthy tissue justifies switching back to 1. A major difficulty is of course the scarce availability of biological markers, which critically depends on each particular cancer. For instance, in prostate cancer, a regrowth of the plasmatic level of PSA, routinely available to clinical measures for quite a long time, after some stagnation time under treatment may indicate the emergence of resis- tance. In the same way, for colorectal cancer, it has been advocated that circulating tumoral DNA detection may be used for clinical manage- ment [54], and the same is true of circulating tu- mour cells [10]; however these techniques are far from being clinical routine. As regards damages to healthy tissue, they are numerous (e.g., for 5-FU and other cytotoxic drugs, classical hand- foot syndrome, mouth sores, neutropenia, that often lead to treatment interruption), depending on each molecule and most of all on the evalu- ation of their severity by the oncologist in the clinic, given the health status of the patient un- der treatment (in the case of laboratory animals, weight loss is a common indicator of toxicity). Taking advantage of the models introduced in Section III, there are several directions for analysing such types of optimal control but in a slighly different or generalised setting. This would both test the robustness of the strategy presented above, and possibly lead to alternative ones depending on the context. The addition of an advection term would help modulate the speed at which emergence and resensitisation occur. Modelling how such terms would depend (or not) on a given drug is already an issue. However, it is likely that the addition of such a term will not jeopardise the numerical computation of optimal controls with direct methods refined with homotopies. Of course, considering higher-dimensional phe- notype variables or adding a space variable will inevitably lead to an explosion in complexity. This is why the numerical optimal control of phenotype-structured PDEs will benefit from state-of-the-art methods in that field, which in turn highly rests on the quality of optimisa- tion solvers. In other words, expert numerical methods will undoubtedly be at the core of any attempt at solving these complex infinite- Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 15 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... dimensional problems. The theoretical aspects are more exploratory. Even for the integro-differential system of [81], the optimal control had to be solved in a re- stricted class of controls, and a complete under- standing of the interplay between concentration phenomena and optimal control is yet to emerge. With few tools available for the control and op- timal control of non-local PDEs (an active area of research), a theoretical analysis of optimal controls of PDE structured models is at this stage a real hurdle. V. Open and challenging questions A. Conflicting phenotypes and multicellularity The question of emergence of drug-resistant clones in a possibly totally genetically homoge- neous cancer cell population (as is likely the case of the observations reported in [76], [77], [88]) under environmental pressure, here drug expo- sure, is related to the emergence of multicellu- larity in unicellular organisms. This question has been the object of many studies by evolutionary biologists [1], [63], [65], [66], [67], and they hypothesise that, confronted with a challeng- ing, possibly deadly, environmental pressure, an already existing, without specialisation, multi- cellular aggregate (this had to occur after the beginning of massive oxygenation of the ocean and atmosphere, about one billion years ago, as, to stick together, cells need some glue of collagen family, which is synthesised only in the presence of free oxygen [94], [95]) had to specialise to survive. The proposed paradigmatic scenario is in [63], [65] the conflict between proliferation - or fecundity, with adhesivity, to maintain against predators a colony of replicating cells on a good environmental trophic niche - and motility - to make the aggregate able to change its location, to leave for a more favourable one when re- sources are become scarce or when predators are threatening. The solution of such conflict is found in specialisation in two phenotypes, later to be refined, likely by bifurcations in more than two if the environmental pressure is diversified. This question has been tackled by Yannick Viossat together with Richard Michod in a simple setting [67], from which one finds that according to the convexity or concavity of a level set on which an optimum of fitness is to be found, there may be coexistence of two phenotypes or predominance of a single one. Such a situation is encountered in an adaptive dynamics framework in [70] (see Section III-B) for the possible invasion, or coexistence with healthy cells, of a leukaemic cell clone. However, how to model such specialisation and cooper- ativity in general, or in the particular case of viability vs. fecundity for cancer cells under drug pressure, still escapes our efforts. B. Coherence in an organism and its control The question of coherence - and of within- and between-tissue cohesion - of a whole mul- ticellular organism with so many diverse and specialised subpopulations is seldom posed, and except in [78], it is generally ignored. Matej Plankar and co-authors ask the main question that is so often dodged when speaking of cancer as a developmental disorder: ‘What exactly is disorganised that was previously organised?’ and they propose that the biophysical base of such coherence resides in the coherence, in the phys- ical sense, of oscillatory signals between cells that might be of electromagnetic or quantum nature, transporting energy and information, that could originate from, or be transmitted by, microtubules working like antennas, to - likely too vaguely and unfaithfully - sum up their hypotheses. How are such signals synchronised? Is there a forcing signal originating from an organisational centre, or is it based on some sort of multilevel system of phase-locked loops? A possible candidate for such an organising system is the circadian system, that is made of circadian clocks [84], [87], existing in all nucle- ated cells (and even, so it seems, in red blood cells by different mechanisms), and that consist of oscillators based on gene networks that exist in all, at least, animal cells (but have also been individuated in some plants). Such oscillators Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 16 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... have firstly been evidenced in fruitflies [51], and later in all mammals [105] where they have been searched for. They have the general property to be daily reset by the sun (or by routine social activities when the sun does not shine its rays), and they date back to a very ancient past of our planet. There exists a central control centre, the circadian pacemaker located in the suprachias- matic nuclei of the hypothalamus in mammals, that receives synchronising electric signals from external light via the retinohypothalamic tract and physiologically sends synchronising mes- sages to all peripheral cells via hormones and the autonomic nervous system. The activity of the central circadian pacemaker, that is in particular reflected by body temperature oscillations and by oscillations of corticosteroids in the surrenal gland, is known to be disrupted in cancer [85], and this all the more so as cancer is more advanced. Although the authors of [78] do not mention this synchronising system, it is coherent with their view. Is the circadian the synchronis- ing system? or is it a dubbing system, under the dependence of an electromagnetic or quan- tum signalling system advocated in [78]? More hidden than the circadian system, could some organising coded plan, progressively established together with the immune system in the devel- opment of multicellularity, be the glue and con- trol on which all Metazoan between- and within- tissue coherence relies? In other words, could there exist a set of genes, already present in early Metazoa, likely related to epigenetic control, that would on the one hand define, in a MHC (major histocompatibility complex)-like way, in its fixed part a species and, at the individual level, an individual within a species, and on the other hand a variable part within a species that would give rise to the different cell phenotypes that make a coherent multicellular organism (in limited number, 200 to 400 cell types or so, from enterocytes to neurons in a Human). Would this be the case, then one could imagine that tumours - as Metazoa 1.0, according to the atavistic hypothesis of cancer - might have developed a sort of primitive, failed, immune re- sponse system whose main failure and difference with respect to the host normal immune system would be a strong tolerance to plasticity, i.e., to lack of belonging to a well-differentiated cell class. In the metaphoric Waddingtonian view, this would imply an ablated, flattened epigenetic landscape, with plenty of room for dediffer- entiation and transdifferentiation between cell phenotypes. Could such flattened Waddington’s epigenetic barriers in return be interpreted in terms of undecided, empty, spins borne on cell antigens that should normally, to avoid recogni- tion as foe by antigen presenting cells, be coded as either 1 (differentiated) or 0 (open to further differentiation in a well-defined cell fate), but not blank? Some support to these speculations may be found in a study dedicated to the origin of the Metazoan immune system [68]. C. Intra-tumour cooperativity, plasticity, bet hedging If tumours, as Metazoa 1.0, have developed some internal cooperativity that makes them able to survive as a whole to cytotoxic stress and to friend-or-foe recognition by the immune system, what does such cooperativity consist of? Experimental evidence exists that such cooper- ativity exists [21], [79], [92] and is necessary for a tumour to thrive, while some studies focus on evidencing tumour genetic or only pheno- typic heterogeneity [61], [62] without necessarily proposing a rationale for such heterogeneity. However, as pointed by Mark Vincent, “Hetero- geneity, even though it might in some superficial way, ‘explain’ differential drug sensitivity, is not in itself an explanation of cancer; rather, it is the heterogeneity itself that requires explana- tion.” [104]. Indeed, focusing only on drug resis- tance, we might satisfy ourselves with the plain observation of tumour heterogeneity to explain the variety of drug resistance mechanisms and why it is so uneasy to eradicate cancer. But try- ing to understand its determinants is more diffi- cult. Leaving aside the obvious fact that spatial isolation of cells inside a tumour may lead under Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 17 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... different forms of environmental pressure to dif- ferent (phenotypically or genetically) clones that may have nothing to do with cooperativity, we may wonder why different phenotypes may be found in the same spatial niche. To begin with, does cooperativity exist with a fixed repartition of phenotypes inducing some division of labour in a tumour, or is it not a transient state that is observed only when a cancer cell population is put at stake under cellular stress?... Or in artificial lab conditions [21], [79], [92]? Can we consider that no actual cooperativity, in the sense of division of labour in an integrated structure, exists within a cancer cell population, all the more so as cell differentiations are im- paired in cancer, but that plasticity of cancer cells (and not only of cancer populations) is so high - within a preestablished Metazoa 1.0 plan, i.e., it is not infinite, but takes advan- tage of many, but finitely many, stress response mechanisms inscribed in their genome and easily reactivatable - that tumours can react to deadly insults by different resistance mechanisms, the simplest one being enhanced proliferation out of control, making them winners in all (known to their genome) cases? The sole idea of cancer cooperativity should be examined with care, if one admits that cancer cells are fundamentally cheaters, or otherwise said, defectors in the evolutionary game of multicellularity [2]. How- ever, primitive Metazoa, such as sponges [68] or algae show some cooperativity, in particular as regards immunity to invasion by pathogens. Have successful tumours regressed in evolution at an earlier stage than sponges? Likely yes, as multicellularity in sponges is well controlled. Following the theme of cancer cell plastic- ity, an interesting notion has recently emerged, the so-called bet hedging fail-safe strategy of cancer cell populations [9], [41]. According to this hypothesis, cancer cells - or cancer cell populations - are so plastic that they can adapt their phenotypes to sustain different insults in- volving critical cell stress by developing differ- ent adapted subpopulations. It has also been observed that some cancer cells may express very ancient genes (so-called ‘cold genes’, i.e., that are conserved throughout evolution, being protected from evolution due to their essential role in facing unpredictable, but already met in a remote past of evolution, deadly insults) in case of exposure to chemotherapies [107]. One could speculate that some sentinel cells, expressing these ‘cold genes’, might send various resistance messages to other cells, or that they could them- selves, being extremely plastic, differentiate into diverse categories of cell subpopulations, each one developing one of the resistance mechanisms elaborated in the course of evolution from a protozoan state, and then sheer darwinian se- lection would prevail. Whether cells themselves are plastic and can adapt in a sort of Lamar- ckian (necessarily epigenetic) way or only cell populations are plastic, constituted by preexist- ing (prior to any insult) genetically well-defined subpopulations is not easy to decide, and in [13], both scenarios were challenged by a reaction- diffusion-advection model (see Section III-D). However, the very fact that cell differentiations are always - to some extent - impaired in can- cer cells, the fact that inhibitors of epigenetic enzymes have been shown in some cases to an- nihilate drug resistance in very aggressive forms of cancer (NSCLC cells in culture in [88]) induce us to propose plasticity as a distinctive and continuous trait of cancer cells. With respect to the class of cell-functional models proposed in Section III-D, i.e., structured in the conflicting continuous traits named viability (x, potential of survival in extreme conditions, opposed to apop- tosis) and fecundity (y, proliferation potential), one could then add a plasticity trait (θ, opposed to differentiation), characterising, together with the first two, each cell in its relevant variability in a heterogeneous cancer cell population (which would not give an explanation of such hetero- geneity, but might help understand its evolution under cellular stress). A general class of cell population adaptive dynamics models that would be structured in Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 18 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... (x,y,θ) could then, following an idea popu- larised in [6] for the so-called cane toad equation (that describes the invasion of cane toads in Australia by using an equation than cannot be of the classical reaction-diffusion type yielding travelling waves), be described by nt + ∇·{V (x,y,θ,D) n} = α(θ)nxx + β(θ)nyy + γ(θ)nθθ +n { r(x,y,θ) −µ(x,y,θ,D) − ρ(t) C(x,y) } , where r(x,y,θ) is the intrinsic (i.e., in the absence of any limitation) growth rate of the population, µ is an added death term due to the drug dose D, condition C(x,y) ≤ K represents an environmental constraint, V an optional ad- vection function standing for abrupt modifica- tions of the environment (such as the major cell stress-inducing delivery of high doses of drugs as in [13]) and ρ(t) = ∫ [0,1]3 n(x,y,θ,t) dxdy dθ is the total cell population at time t, put as usual in Lotka-Volterra settings in a logistic position to represent competition, e.g., for nu- trients, hence growth limitation, within the cell population. How such class of models might lead to repre- sent the emergence of different cell subpopula- tions under environmental pressure is still work underway. VI. Conclusion In this review of models of adaptive dynam- ics dedicated to represent, analyse and control drug-induced drug resistance in cancer, we have firstly proposed a brief description of the biolog- ical background of cancer evolution, describing in particular the atavistic hypothesis of cancer, which to our meaning illuminates the scenery of the cancer disease, with more and more observa- tion facts to support it. Then (neglecting clas- sical compartmental ODE models that cannot claim to represent adaptive phenomena), we pre- sented different cell population dynamics mod- els that belong to the mathematical category of adaptive dynamics, i.e., integro-differential or partial differential equations structured by continuous traits describing the heterogeneity of cancer cell populations and their evolution under drug exposure. In a third part, we showed how optimal control methods can be applied to such equations of adaptive dynamics and used to design theoretical optimal therapeutic control strategies. Such strategies, even though methods for the identification of the parameters and functions of the models still remain to be found, are amenable to predict qualitative behaviour of cancer cell populations under optimised time- scheduled drug exposure. Finally, we presented some challenging questions, addressed to evolu- tionary biologists and ecologists of cancer, to on- cologists, and to mathematicians to accurately represent, analyse and control the behaviour of cancer cell populations. Acknowledgments The authors are gratefully indebted to- wards their colleagues Luis Almeida, Rebecca Chisholm, Tommaso Lorenzi, Alexander Lorz, Benoît Perthame, and Emmanuel Trélat, co- authors with them of the mathematical pa- pers on structured cell population models of evolution towards drug-induced drug resistance, whose main results have been presented in this review. C.P. acknowledges support from the Swedish Foundation of Strategic Research grant AM13-004. References [1] A. Aktipis, et al., Life history trade-offs in cancer evolution, Nature Rev. Cancer (2013), 13:883-892. [2] A. Aktipis, et al., Cancer across the tree of life: cooperation and cheating in multicellularity, Phil. Trans. R. Soc. B (2015) 370:20140219. [3] L. Almeida, et al., Evolution of cancer cell pop- ulations under cytotoxic therapy and treatment optimisation: insight from a phenotype-structured model, ESAIM Math. Model. Numer. Anal. (2018), accepted. Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 19 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... [4] D. Basanta, et al., Investigating prostate cancer tumour-stroma interactions: clinical and biological insights from an evolutionary game. British Journal of Cancer (2011), 106(1):174-181. [5] C. Basdevant, et al., Optimisation of time- scheduled regimen for anti-cancer drug infu- sion, ESAIM Math. Model. Numer. Anal. (2005), 39(6):1069-1086. [6] E. Bouin, et al., Invasion fronts with variable motil- ity: Phenotype selection, spatial sorting and wave acceleration, C. R. Acad. Sci. Paris, Ser. I (2012), 761-766. [7] T. Boveri, The origins of malignant tumors, Williams & Wilkins, Baltimore (1929). [8] N.F. Britton, Essential Mathematical Biology, Springer, 2003. [9] B. Brutovsky & D. Horvath, Structure of Intra- tumor Heterogeneity: Is Cancer Hedging Its Bets? arXiv (2013),1307.0607. [10] C. Burz, et al. Circulating tumor cells in clinical research and monitoring patients with colorectal cancer, Oncotarget (2018), 9(36):24561-24571. [11] K.J. Bussey, et al., Ancestral gene regulatory net- works drive cancer, PNAS (2017), 114(24):6160- 6162. [12] C. Carrère, Optimization of an in vitro chemother- apy to avoid resistant tumours, J. Theor. Biol. (2017), 413:24-33. [13] R.H. Chisholm, et al., Emergence of drug tol- erance in cancer cell populations: an evolution- ary outcome of selection, non-genetic instability and stress-induced adaptation. Cancer Research, 75(6):930-939, 2015. [14] R.H. Chisholm, et al., Cell population heterogene- ity and evolution towards drug resistance in can- cer: Biological and mathematical assessment, theo- retical treatment optimisation, Biochem. Biophys. Acta (2016), 1860:2627-2645. [15] R.H. Chisholm, et al., Effects of an advection term in nonlocal Lotka-Volterra equations. Comm. Math. Sciences, 14:1181-1188, 2016. [16] H. Cho & D. Levy, Modeling continuous levels of resistance to multidrug therapy in cancer. Appl. Math. Model. (2018), 64:733-751. [17] H. Cho & D. Levy, Modeling the chemotherapy- induced selection of drug-resistant traits during tumor growth. J. Theor. Biol. (2018), 436 (7):120- 134. [18] P. Cirri & P. Chiarugi, Cancer associated fibrob- lasts: the dark side of the coin, Am. J. Cancer Research (2011), 1(4):482-497. [19] L.H. Cisneros, et al., Ancient genes establish stress- induced mutation as a hallmark of cancer. PLoS One (2017) 12(4):e0176258. [20] J. Clairambault & O. Fercoq, Physiologically struc- tured cell population dynamic models with ap- plications to combined drug delivery optimisation in oncology. Mathematical Modelling of Natural Phenomena (2016), 11(6):45-70. [21] A.S. Cleary, et al., Tumour cell heterogeneity main- tained by cooperating subclones in Wnt-driven mammary cancers, Nature Lett. (2014), 508:113- 117. [22] M.I.S. Costa, et al., Optimal chemical control of populations developing drug resistance, Mathemat- ical Medicine and Biology (1992), 9(3):215-226. [23] P.C.W. Davies & C.H. Lineweaver, Cancer tumors as metazoa 1.0: tapping genes of ancient ancestors, Phys. Biol. (2011), 8(1):015001. [24] L. Desvillettes et al., On selection dynamics for con- tinuous structured populations, Commun. Math. Sci. (2008), 6(3):729-747. [25] L. Ding, et al., Clonal evolution in relapsed acute myeloid leukaemia revealed by whole genome se- quencing, Nature (2012) 481:506-510. [26] B. Dirat, et al., Cancer-associated adipocytes ex- hibit an activated phenotype and contribute to breast cancer invasion, Cancer Research (2011), 71(7):2455-2465. [27] T. Dobzhansky, Biology, molecular and organismic, Am Zool. (1964), 4:443-452. [28] T. Domazet-Lošo & D. Tautz, An ancient evolu- tionary origin of genes associated with human ge- netic diseases, Mol. Biol. Evol. (2008), 25(12):2699- 2707. [29] T. Domazet-Lošo & D. Tautz, Phylostratigraphic tracking of cancer genes suggests a link to the emer- gence of multicellularity in metazoa, BMC Biol. (2010), 8(1):66. [30] H. Easwaran, et al., Cancer Epigenetics: Tumor Heterogeneity, Plasticity of Stem-like States, and Drug Resistance, Molecular Cell (2014), 54:716- 727. [31] A.P. Feinberg, et al., Epigenetic modulators, modi- fiers and mediators in cancer aetiology and progres- sion, Nature Rev. Genet. (2016), 17(5):284-99. [32] R.A. Gatenby, et al., Adaptive therapy. Cancer Research (2009), 69(11):4894-4903. [33] M. Gerlinger, et al., Intratumor heterogeneity and branched revealed by multiregion sequencing. N Engl J Med (2012), 366(10):883-892. [34] M. Gerlinger & C. Swanton, How darwinian models inform therapeutic failure initiated by clonal het- erogeneity in cancer medicine. Br J Cancer (2010), 103(8):1139-1143. [35] R.J. Gillies, et al., Evolutionary dynamics of car- cinogenesis and why targeted therapy does not work, Nature Rev. Cancer (2012), 12(7):487-493. [36] A. Goldman, et al., Integrating Biological and Mathematical Models to Explain and Overcome Drug Resistance in Cancer, Part 1: Biological Facts and Studies in Drug Resistance, Current Stem Cell Reports (2017), 3:253-259. [37] A. Goldman, et al., Integrating Biological and Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 20 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... Mathematical Models to Explain and Overcome Drug Resistance in Cancer, Part 2: From Theoreti- cal Biology to Mathematical Models, Current Stem Cell Reports (2017), 3:260-268. [38] M.M. Gottesman. Mechanisms of cancer drug re- sistance. Annu Rev Med (2002), 53:615-627. [39] M. Greaves, Cancer stem cells: back to Darwin? Semin Cancer Biol (2010), 20(2):65-70. [40] M. Greaves & C.C. Maley. Clonal evolution in cancer, Nature (2012), 481(7381):306-313. [41] C.A. Gravenmier, et al., Adaptation to Stochastic Temporal Variations in Intratumoral Blood Flow: The Warburg Effect as a Bet Hedging Strategy, Bull Math Biol (2017), 80(5):954-970. [42] G.D. Guler, et al., Repression of Stress-Induced LINE-1 Expression Protects Cancer Cell Subpop- ulations from Lethal Drug Exposure, Cancer Cell (2017), 32:221-237. [43] P. Hirsch, et al., Genetic hierarchy and temporal variegation in the clonal history of acute myeloid leukaemia, Nature Comm. (2016) 7:12475. [44] S. Huang, et al., Bifurcation dynamics in lineage- commitment in bipotent progenitor cells, Dev Biol, (2007), 305(2):695-713. [45] S. Huang, On the intrinsic inevitability of cancer: From foetal to fatal attraction, Sem. Canc. Biol. (2011), 21:183-199. [46] S. Huang, Genetic and non-genetic instability in tu- mor progression: link between the fitness landscape and the epigenetic landscape of cancer cells, Canc. Metastasis Rev. (2013), 32:423-448. [47] L. Israel, Tumour progression: random mutations or an integrated survival response to cellular stress conserved from unicellular organisms? J. Theor. Biol. (1996), 178(4):375-380. [48] P.E. Jabin & G. Raoul, On selection dynamics for competitive interactions, J. Math. Biol. (2011), 63:493-551. [49] F. Jacob, Evolution and tinkering, Science (1977), 196 (4295):1161-1166. [50] M. Kimmel & A. Świerniak, Control Theory Ap- proach to Cancer Chemotherapy: Benefiting from Phase Dependence and Overcoming Drug Resis- tance, in Tutorials in Mathematical Biosciences III, Avner Friedman ed., Springer, 2006. [51] R.K. Konopka & S. Benzer, Clock mutants of drosophila melanogaster, Proc Natl Acad Sci USA (1971), 68:2112-2116. [52] O. Lavi, et al., The role of cell density and in- tratumoral heterogeneity in multidrug resistance, Cancer Res (2013), 73(24):7168-7175. [53] U. Łędżewicz, et al., On the role of tumor hetero- geneity for optimal cancer chemotherapy, Networks & Heterogeneous Media (2019), 14(1):131-147. [54] H. Li, et al., Circulating tumor DNA detection: A potential tool for colorectal cancer management, Oncology Letters (2019), 7(2):1409-1416. [55] C.H. Lineweaver, et al., Targeting cancer’s weak- nesses (not its strengths): therapeutic strategies suggested by the atavistic model, Bioessays (2014), 36(9):827-835. [56] T. Lorenzi, et al., Dissecting the dynamics of epige- netic changes in phenotype-structured populations exposed to fluctuating environments. Journal of Theoretical Biology, 386:166-176, 2015. [57] T. Lorenzi, et al., Tracking the evolution of can- cer cell populations through the mathematical lens of phenotype-structured equations. Biology Direct, 11:43, 2016. [58] A. Lorz, et al., Populational adaptive evolution, chemotherapeutic resistance and multiple anti- cancer therapies, ESAIM Math. Model. Numer. Anal. (2013), 47(2):377-399. [59] A. Lorz, et al., Modeling the effects of space struc- ture and combination therapies on phenotypic het- erogeneity and drug resistance in solid tumors, Bull. Math. Biol. (2014), 77(1):1-22. [60] Y. Manabe, et al., Mature adipocytes, but not preadipocytes, promote the growth of breast carci- noma cells in collagen gel matrix culture through cancer-stromal cell interactions, The Journal of pathology (2003), 201(2):221-228. [61] A. Marusyk, et al., Intra-tumour heterogeneity: a looking glass for cancer?, Nature Rev. Cancer (2012), 12: 323-334. [62] A. Marusyk, et al., Non-cell-autonomous driving of tumour growth supports sub-clonal heterogeneity, Nature (2014), 514:54-58. [63] J. Maynard Smith & Eös Szathmáry. The Major Transitions in Evolution. Oxford University Press, 1997. [64] K.D. McCullough, et al., Plasticity of the neoplastic phenotype in vivo is regulated by epigenetic factors. Proc Natl Acad Sci USA (1998), 95(26):15333- 15338. [65] R.E. Michod & D. Roze, Cooperation and conflict in the evolution of multicellularity, Heredity (2001), 36:1-7. [66] R.E. Michod & D. Roze, Cooperation and conflict during evolutionary transitions in individuality, J. Evol. Biol. (2006), 19:1406-1409. [67] R.E. Michod, et al., Life-history evolution and the origin of multicellularity, J. Theor. Biol. (2006), 239:257-272. [68] W.E.G. Müller & I.M. Müller, Origin of the Metazoan Immune System: Identification of the Molecules and Their Functions in Sponges, Inte- grative Comparative Biology (2003), 43:281-292. [69] J.D. Murray, Mathematical biology, 2 volumes, Springer, 2002, 2003. [70] T.N. Nguyen, et al., Adaptive dynamics of hematopoietic stem cells and their support- ing stroma: A model and mathematical anal- Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 21 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... ysis, Submitted (2019), preprint available at https://hal.archives-ouvertes.fr/hal-01963820. [71] P.C. Nowell, The clonal evolution of tumor cell populations. Science (1976), 194(4260):23-28. [72] A. Olivier & C. Pouchol, Combination of direct methods and homotopy in numerical optimal con- trol: application to the optimization of chemother- apy in cancer, Journal of Optimization Theory and Applications (2018), online Dec. 18, 2018, 25 pages. [73] E. Pasquier, et al., Metronomic chemotherapy: new rationale for new directions, Nature Rev. Clin. Oncol. (2010), 7(8):455-465. [74] B. Perthame, Transport equations in biology, Springer, 2007. [75] B. Perthame, Parabolic equations in biology, Springer, 2015. [76] A.O. Pisco, et al., Non-darwinian dynamics in therapy-induced cancer drug resistance, Nat. Com- mun (2013), 4:2467. [77] A.O. Pisco & S. Huang, Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: ‘What does not kill me strengthens me’, Br J Cancer (2015), 112(11):1725-1732. [78] M. Plankar, et al., On the origin of cancer: Can we ignore coherence? Progress in Biophysics and Molecular Biology (2011), 106(2):380-390. [79] K. Polyak & A. Marusyk, Clonal cooperation, Na- ture (2014), 508:52-53. [80] C. Pouchol, Modelling interactions between tu- mour cells and supporting adipocytes in breast cancer (2015), Internship report available at https://hal.inria.fr/hal-01252122. [81] C. Pouchol, et al., Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy, J. Math. Pures Appl. (2018), 116:268-308. [82] C. Pouchol & E. Trélat, Global stability with selec- tion in integro-differential Lotka-Volterra systems modelling trait-structured populations, Journal of Biological Dynamics, 12:1, 872-893, 2018. [83] C. Pouchol, Analysis, control and optimisation of PDEs, application to the biology and ther- apy of cancer, PhD thesis, Sorbonne Université (June 2018), available at https://hal.inria.fr/tel- 01889253. [84] S.M. Reppert & D.M. Weaver, Coordination of cir- cadian timing in mammals, Nature (2002), 418:935- 941. [85] T. Rich, et al., Elevated serum cytokines correlated with altered behavior, serum cortisol rhythm, and dampened 24-hour rest-activity pattern in patients with metastatic colorectal cancer, Clinical Cancer Research (2005), 11:1757-1764. [86] H. Schättler & U. Łędżewicz, Optimal Control for Mathematical Models of Cancer Therapies, Springer, 2015. [87] U. Schibler & P. Sassone-Corsi, A web of circadian pacemakers, Cell (2002), 111:919-922. [88] S.V. Sharma, et al., A chromatin-mediated re- versible drug-tolerant state in cancer cell subpopu- lations, Cell (2010), 141(1):69-80. [89] L. Sigalotti,et al., Epigenetic drugs as immunomod- ulators for combination therapies in solid tumors, Pharmacology & Therapeutics (2014), 142:339-350. [90] A.M. Soto & C. Sonnenschein, The somatic muta- tion theory of cancer: growing problems with the paradigm? BioEssays (2004), 26(10):1097-1107. [91] G.W. Swan, Role of Optimal Control Theory in Cancer Chemotherapy, Mathematical Biosciences (1990), 101:237-284. [92] D.P. Tabassum & K. Polyak, Tumorigenesis: it takes a village, Nature Rev. Cancer (2015), 15:473- 483. [93] G. Tonini, et al., Rechallenge therapy and treat- ment holiday: different strategies in management of metastatic colorectal cancer, Journal of experi- mental & clinical cancer research (2013), 32(1):1. [94] K.M. Towe, Oxygen-collagen priority and the early metazoan fossil record, Proc Natl Acad Sci USA (1970), 65(4):781-788. [95] K.M. Towe, Biochemical keys to the emergence of complex life, p. 297-305 in J. Billingham ed., Life in the Universe, MIT Press, Cambridge, MA, 1981. [96] E. Trélat, Optimal control and applications to aerospace: some results and challenges, Journal of Optimization Theory and Applications (2012), 154(3):713-758. [97] A.S. Trigos, et al., Altered interactions between unicellular and multicellular genes drive hallmarks of transformation in a diverse range of solid tumors, PNAS (2017), 114 (24):6406-6411. [98] A.S. Trigos, et al., How the evolution of multicellu- larity set the stage for cancer, Br. J. Cancer (2018), 118:145-152. [99] J.E. Trosko, Mechanisms of tumor promotion: pos- sible role of inhibited intercellular communication, Eur. J. Cancer Clin. Oncol. (1987), 23(6):599-601. [100] J.E. Trosko, Gap Junctional Intercellular Commu- nication as a Biological “Rosetta Stone” in Under- standing, in a Systems Biological Manner, Stem Cell Behavior, Mechanisms of Epigenetic Toxicol- ogy, Chemoprevention and Chemotherapy, J. Mem- brane Biol. (2007), 218:93-100. [101] J.E. Trosko, A conceptual integration of extra-, intra- and gap junctional-intercellular communica- tion in the evolution of multi-cellularity and stem cells: how disrupted cell-cell communication during development can affect diseases later in life, Int. J. Stem Cell Res. Ther. (2016), 3(1):021. [102] B. Ujvari, B. Roche, F. Thomas, Eds., Ecology and evolution of cancer, Elsevier Academic Press, 2017. [103] M.D. Vincent, Cancer: a de-repression of a default Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 22 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 J. Clairambault, C. Pouchol, A survey of adaptive cell population dynamics models of emergence ... survival program common to all cells?: a life-history perspective on the nature of cancer, Bioessays (2011), 34(1):72-82. [104] M.D. Vincent, Cancer beyond speciation, Ad- vances in cancer research (2011), 112:283-350. [105] M.H. Vitaterna, et al., Mutagenesis and mapping of a mouse gene, clock, essential for circadian be- havior, Science (1994), 264:719-725. [106] C.H. Waddington, The strategies of the genes, George Allen & Unwin, London (1957). [107] A. Wu, et al., Ancient hot and cold genes and chemotherapy resistance emergence, Proc. Nat. Acad. Sci. USA, (2015), 112:10467-10472. [108] J.X. Zhou, et al., Phylostratigraphic analysis of tumor and developmental transcriptomes reveals relationship between oncogenesis, phylogenesis and ontogenesis, Converg. Sci. Phys. Oncol. (2018), 4:025002. Biomath 8 (2019), 1905147, http://dx.doi.org/10.11145/j.biomath.2019.05.147 Page 23 of 23 http://dx.doi.org/10.11145/j.biomath.2019.05.147 Introduction: Motivation from and focus on drug resistance in cancer Biological background The many facets of drug resistance in cancer Ecology, evolution and cancer in cell populations The atavistic theory of cancer Models of adaptive dynamics Models structured in resistance phenotype Modelling mutualistic tumour-stroma interactions Models structured in phenotype and space Models structured in cell-functional variables Optimisation and optimal control ODE models and their optimal control in cancer Optimal control of phenotype-structured PDE models Future prospects in optimal control Open and challenging questions Conflicting phenotypes and multicellularity Coherence in an organism and its control Intra-tumour cooperativity, plasticity, bet hedging Conclusion References