Mathematical Modelling of Cancer Growth and Drug Treatments: Taking Into Account Cell Population Heterogeneity and Plasticity

Mathematical models of cancer growth and evolution of cancer cell characteristics, aka traits or phenotypes, together with optimisation and optimal control methods to contain them, in the framework of adaptive dynamics of cell populations, are presented. They take into account the heterogeneity of cancer cell populations, i.e., their biological variability, and their intrinsic plasticity, i.e., their nongenetic instability that allows them to quickly adapt to changing environments. The presented vision of the cancer disease, which is specific to multicellular organisms, relies on a relatively novel vision, consistent with a billion-year evolutionary perspective. Based on recent contributions from philosophy of cancer, these mathematical models aim at designing theoretical therapeutic strategies to simultaneously contain tumour progression and limit adverse events of drugs to healthy cell populations.


I. INTRODUCTION
"What is cancer?" is a question that must somehow be answered with clear ideas before proposing relevant mathematical models of cancer growth and phenotype evolution, together with methods to optimise its treatments.Recent contributions from philosophy of biology, immunology, and archeology of genomes, are reported here, stressing the importance of a default of cohesion within and between tissues of a same organism -as cancer is a disease of multicellular organisms only -, that is manifested as an impairment of organismic control on both proliferation and differentiation.
Different mathematical models of tumour growth involving cancer cell population heterogeneity and plasticity, and drug treatments of cancer are presented in the present paper.The ideas developed here are not a review of the vast domain of studies dedicated to modelling heterogeneous cell populations.Devoted to advocating and illustrating by a few chosen models the interest of continuous phenotypestructured equations in cancer and its treatment, they are intended to be presented as part of a tutorial session in ECC 2023 on "Modelling and Treatment of Cancer".They concern deterministic models, with a focus set on adaptive dynamics of cell populations, healthy and cancer, and the effects on them of anticancer drugs.Constraints linked to avoiding the main pitfalls of therapeutics in oncology, namely therapy resistance in cancer cells and adverse toxic events in healthy cells, together with an objective function that is most often limited to maximising tumour cell kill, are considered, paving the way for the use of optimisation methods.Optimised clinical strategies for the continuous delivery of anticancer drugs may thus be designed by resorting to optimal control methods.Such optimised strategies should thwart some known pathological strategies used by cancer cell populations to dodge deadly insults induced by drugs, such as the development of drug-resistant cell phenotypes and the so-called 'bet hedging' by which tumours may simultaneously face different modes of insults in abruptly changing environments [5], [23], [46].
From a modelling point of view, the complete reversibility of the cell phenotypes under study when the environment changes is represented by the reversible behaviour of the solutions of the equations in terms of cell-structuring traits (aka phenotypes).It is consistent with the fact that we study in an isogenic context cell differentiations and dedifferentiations, i.e., reversal of differentiations within a given cell lineage toward a less mature phenotype [54], that rely on epimutations, i.e., epigenetic changes in the expression of genes, and not on gene mutations [24].With or without interactions between cell populations, the equations are of the nonlocal Lotka-Volterra type, structured by adaptive continuous cell phenotypes, with optional diffusion and advection terms (two examples are given below).Adaptation, i.e., phenotype changes under environmental pressure, e.g., due to competing cell populations, drugs or other therapeutic means, is thus considered at a cell population status at which the situation of a cancer cell population is reversible (i.e., not after mutations), in particular for the occurrence of druginduced drug resistance.This allows to theoretically justify the clinical notions of drug holiday, i.e., stopping delivering a drug for a defined time interval, and re-challenge, i.e., giving the same drug -and not another one by a wrong consideration that the first one has definitely failed or will fail -after a drug holiday, in cancer [56].

II. WHERE BIOLOGY MEETS PHILOSOPHY A. Heterogeneity and plasticity of cancer cell populations
Heterogeneity of cancer cell populations [37] has often been mentioned as a major obstacle to the efficacy of drug treatments, in particular of so-called 'targeted therapies' [22], as treatments that may be very effective on a dominant genetic clone of a tumour cease to be so when another clone emerges [15], be it in malignant haemopathies [15] (aka liquid tumours) or in solid tumours [21].This certainly concerns genetic heterogeneity, i.e., situations in which spontaneous mutations have appeared due to the genetic instability of tumour cells, but isogenic heterogeneity resulting from nongenetic instabilities in gene expression due to the lability of epigenetic enzyme activity [53] is as much important [19], [20], [38], and it occurs before mutations, that remain rare events.The plasticity of cancer cells with respect to the expression of genes [12], [19], [38], [60], as all differentiations, that make an isogenic multicellular organism so varied in cell types (20 for sponges [41], [42], about 200 for human beings) is due to epigenetic factors.These are strictly controlled at the organismic level in physiology, and out of control in cancer [19], [20].Such plasticity is characterised by the velocity of adaptation of cell phenotypes, such as fecundity, viability or motility, to changes in the cell, tissue and whole organism ecosystem.

B. Philosophy of cancer
A recently autonomised field of knowledge, the so-called philosophy of cancer [4], [12], [55] considers cancer as an impairment of organismic control on both proliferation and differentiations [4].It should be stressed that cancer is a disease that can be met only in the animal kingdom, constituted of multicellular organisms, that, different from plants (plants may also develop tumours, but they remain localised in capsules, never threatening the existence of the organism), are endowed with other properties such as motility and heterotrophy.Multicellular animals, also called Metazoa, emerge from a unique cell, the zygote (fecundated egg) and are thus completely isogenic.The same genome thus exists in all cells of a given animal organism, and the extreme diversity of cell types that constitute an animal is only due to successive differentiations, by selective expression or repression of genes at each differentiation step.Such expressions or repressions are normally under strict epigenetic control, i.e., under the control of epigenetic enzymes that may become impaired in cancer.Furthermore, in cancer, these epigenetic mechanisms are easily manipulated in the plastic cancer cell populations in response to changes of the tumour microenvironment (such as due to the introduction of life-threatening drugs) to develop (reversible) resistant cell phenotypes.The question of understanding what makes a multicellular individual coherent, contrary to a tumour, may not always be summed up by its zygotic origin, and it is the object of the philosophy of immunology [50], [51], as one of the components that define a multicellular individual is the necessary existence of an immune system [41], [42], absent in tumours, which makes possible the success of immunotherapies.

C. The atavistic theory of cancer
In a completely independent way from the above mentioned works, the atavistic theory of cancer was speculated and popularised from 2011 on in a series of papers [13], [58], [31], relying or partly confirmed in previous [25], or following papers examining phylostratigraphic trees [17], [18], [32], [33] or archeogenetic data [6], [10], [57].It proposes that cancer is a reversal in evolution from a coherent multicellular organism towards an incoherent, nonviable state closer to unicellularity, but representing an autonomous, coarse and unsuccessful try towards multicellularity from partly organised collections of cells that escape organismic control, a state which these authors name Metazoa 1.0.The interest of considering this atavistic theory is that all cancers are constituted of cells in a same multicellular organism, that have kept in the memory of their common genome defence mechanisms against external aggressions from the environment, that can be unleashed and recruited due to the absence of organismic control, that otherwise would thwart them as introducing too loose cohesion of the ensemble, challenging the very existence of the organism [1].Note however that such mechanisms, based on cell de-differentiations and re-differentiations, are at work, strictly controlled, in early development and in adult physiological wound healing.

A. Pitfalls in clinical oncology
In the clinic of cancers, drug treatments meet two main pitfalls: emergence of resistance in cancer cell populations and existence of adverse toxic effects in healthy cell populations.Although most mathematical optimisation works consider one or the other pitfall, it will be shown in the sequel that it is possible to manage them both simultaneously, and not only by following absolute limitations to both population densities, supposed to be taken from clinical considerations, but rather by following with time their respective parts and permanently controlling them.

B. Continuous models of tumour growth
Most models of tumour growth presented until recently were based on ODEs, or PDEs structured in Cartesian space, not many of them are PDEs structured according to phenotypes [35], even fewer according to both phenotypes and space variables [36].The privileged focus of this paper is on structuring tumour growth models in phenotypes [11].

C. Theoretical therapeutic control
Many models of control, including optimal control methods, mimicking therapeutics in the clinic have been proposed, for ODE or PDE models, and useful references on this topic are in particular [26], [30], [52].A particular aspect is control in dormancy [16].In the present article is presented an example of optimal control for a phenotype-structured integro-differential model of the nonlocal Lotka-Volterra type [2], [29], [49].

A. Structured population models and adaptive dynamics
The general idea of structured models is to describe heterogeneity in a population of individuals represented by continuous densities by using continuous traits that are relevant to the question at stake, and that are susceptible to evolve with time by adaptation of such traits to changing environmental conditions or interactions with other populations of individuals.This way of representation originally comes mostly from the seminal paper of McKendrick [39] on public health problems (see in the original paper the equation p. 122) and from theoretical ecology [14], however it has been adapted to cell populations, see, e.g., [45], and a number of following papers.The structuring trait in the population of individuals may be age, size, or physiological traits determining adaptation to the environment such as fecundity, viability [7], [8], [34], motility or plasticity, as will be shown in the sequel.Indeed, as for general and long-term evolution, as regards mutations [40], [43], [59], adaptation to stress induced by the environment, as regards epimutations for short-term non-genetic fast phenotypic changes, should be considered as a major determinant of trait differentiation.

B. Tumours and their micro-environment
Populations of cancer cells, in malignant haemopathies or in solid tumours, are amenable to such representations.Their micro-environment may be constituted of other cell populations, in mutualistic [44] or competitive (in particular tumour-immune, for which many models have been published) interactions.It may also consist of mechanical or chemical constraints, as in particular by the introduction of anticancer drugs, that target different aspects of tumour growth, in particular cell death, intrinsic proliferation or action on the structuring traits.A trait of expression of resistance to a given drug has often been considered [7], [8], [35], [36], [47], [48].

C. Tumour strategies to dodge therapeutic control
As mentioned in the introduction, at least two of these strategies have been identified.In the first one [49], a phenotype of expression of drug resistance structures two cell populations, one tumoral, the other healthy (here healthy only meaning targets of unwanted toxic side effects of the drug), that interact only through a nonlocal logistic term, i.e., an environment integral term I(t), weighted sum of the total mass of the two cell populations, that impinges on the natural death coefficient d(x), representing an assumed competition between all cells for space and nutrients.The model illustrates drug-induced drug resistance, present in cancer, not in healthy cell populations.It will be developed in the next subsection.
The other tumour strategy, that is explored in [3] (without control), is related to tumour bet hedging [5], [23], [46].It hypothesises that in order to adapt to unpredictable changes in the tumour environment, tumour cell populations are amenable to display different cell phenotypes, with differentiated abilities to resist different threats.This is represented in the simplest possible way by phenotypic divergence, i.e., the division of the cell population into two phenotypically distinct subpopulations, each one of them maximising a form of fitness related to two different phenotypes related by a flexible (i.e., not simply y = 1 − x) trade-off condition between them.
The main equation, that describes the evolution of a cell population of density n(t, z) at time t bearing phenotype z, where z = (x, y, θ), with for instance x=viability, y=fecundity, and θ=plasticity, under some environmental pressure represented by an advection term V (t, z), runs as: where (V n − A(θ)∇n) • n = 0 for all z ∈ ∂D (n being the outer normal vector on ∂D) and n(0, z) = n 0 (z) for all z ∈ D = {C(x, y) ≤ K} × [0, 1], D defining a trade-off between traits x and y in [0, 1] (e.g., (x − 1) 2 + (y − 1) 2 ≥ 1), θ free in [0, 1], and the diffusion matrix (where a 11 and a 22 are non decreasing functions of θ) is a representation of how the internal plasticity trait θ impacts the non-genetic instability (due to epigenetic epimutations) of traits x and y, by tuning the diffusion term ∇.{A(θ)∇n}; the advection term where for instance V (t, z) = 10 −3 θ(−y, −x, −(x + y)) in simulations, represents the impact of external evolutionary pressure on the population, by changes in the ecosystem hosting the cell population; the integral ρ(t) = D n(t, z)dz stands for the total mass of cells in the population at time t.
The reader is sent to [3] to check that phenotypic divergence around two different concentrated phenotypes (x 1 , y 1 ) and (x 2 , y 2 ) is shown to actually occur in the model, while the plasticity phenotype θ simultaneously decreases.In this setting where phenotype space is 3-dimensional, the plasticity trait θ tunes the Laplacian term, while the advection function represents the influence of the environment on the adaptable structure variable (3-d phenotype) z = (x, y, θ) present in the cell population with density n(t, z).

D. A phenotype-structured model with optimal control
The model presented in [2] is a consistent illustration of a drug resistance phenotype-structured model in a cancer cell population with optimal control for the optimised delivery of an anticancer chemotherapy acting as an added death term in a reaction-diffusion setting.However, the somewhat richer integro-differential model presented in [49] takes simultaneously into account both emergence of drug-induced drug resistance (to be limited) and adverse toxic effects on a nonresistant healthy cell population.As mentioned in the previous section, the two populations are here linked by mixed population terms I H (t) and I C (t), standing for the influence of a global cell environment in the nonlocal logistic term of the Lotka-Volterra model which is described below.
The phenotype x ∈ [0, 1] of expression of drug resistance (x = 0: total drug sensitivity; x = 1: total drug resistance) structures the two populations.However, the functions r H , d H , µ H for proliferation, death and drug sensitivity, respectively, in healthy cells, are chosen much less sensitive to the drug levels and to the environment variables than their correspondent functions r C , d C , µ C in cancer cells.
(Cancer cell population C) Environment terms involving the two cell populations: Each of these equations for the two cell populations is a transport equation of the nonlocal Lotka-Volterra type, structured by the continuous phenotype x of drug resistance, in which the proliferation term r(x) is tuned by a control term representing the effect of a cytostatic (non-cell killing) drug u 2 , while an added death term µ(x) is tuned by the effect of a cytotoxic (cell killing) drug u 1 .As mentioned earlier, the drug-free death term d(x)I(t) stands for global between-cell competition for space and nutrients.The analysis of the asymptotic behaviour of the system in the absence of drugs (or with constant drug doses) shows concentration of the phenotype x in both populations on a discrete set in [0, 1] [49], which in simulations appears as reduced to a unique point.Then comes the question of optimally controlling this asymptotic behaviour by a combination of the two time-continuous (or L ∞ ) drugs (u 1 , u 2 ).
Optimal control problem OCP: find controls (u 1 , u 2 ) minimising in fixed horizon T the objective function (the last constraint, with, e.g., θ H = 0.6, to limit damage to healthy cells) Then it can be shown that the optimal therapeutic trajectory (u 1 , u 2 ) in large time T > 0 consists of 3 parts: • a long-time part, with constant controls on [0, T 1 ] (where = θ HC ), at the end of which populations have concentrated in phenotype (for T 1 large); • a free arc (no constraint saturating) where controls ) and u 2 = u max 2 .The reader is sent to [49] for proofs, numerical simulations and illustrations.The figure below shows the resulting optimal control (u 1 , u 2 ) for θ HC = 0.4 and θ H = 0.6.

V. DISCUSSION ABOUT STRUCTURED MODELS A. Including targets for control
In all mathematical models, be they phenotype-structured or not, that are designed in the perspective of therapies, Fig. 1.Evolution with time t of the x resistance phenotype distribution curves n H (t, x) and n C (t, x) and solution (u 1 , u 2 ) of the OCP problem for T = 30.In the lower two panels, one can see that the cytotoxic drug u 1 level is maintained positive, but as low as possible, close to zero, for a long time, before delivering its maximal tolerated dose on a short time interval, so as to keep the cell population as sensitive as possible until time T 1 (here T 1 ≈ 20) -which is the drug holiday strategy used in the clinic of cancers [56] -, then delivered at a medium dose.The cytostatic drug u 2 is maintained at a low level, but not close to zero, until some time before and close to T 1 , and then delivered at its maximum tolerated dose.One can see on the two rightmost panels that the cancer cell population ρ C goes to extinction while the healthy cell population ρ H , target of unwanted toxic side effects of the drugs, is preserved above a prescribed level.The central panel shows the evolution of the ratio of drugsensitive cancer cells ρ CS (t) = 1 0 (1 − x)n C (t, x)dx over total cancer cells ρ C (t) = 1 0 n C (t, x)dx.Performed in the modelling language AMPL with expert simulation routines IpOpt.Adapted from [49] with permission.
physiologically based targets of control must be included in the equations.If the control primarily targets cell death (case of most classic chemotherapies), then it must be present in the equations as an added death term.If it targets proliferation, as for, e.g., tyrosine kinase inhibitors (TKIs) that are also antagonists of epidermal growth factor receptors (EGFR antagonists), then it must be present as tuning the intrinsic (that is, defined by setting death terms to zero) instantaneous growth term ∂n n∂t .If it is an epigenetic drug that controls nongenetic instability, then it must address a term tuning a Laplacian ∂ 2 n ∂x 2 or other terms representing nongenetic instability in a reaction-diffusion equation.Controlling the tumour microenvironment might also be thought of as tuning an added advection term ∂n ∂x , the structure variable x being here assumed to be one-dimensional for simplicity of presentation.

B. Taking into account heterogeneity
Various ODE models claim to take account of cell population heterogeneity by compartments representing discrete heterogeneity, such as (totally) drug-sensitive and (totally) drug-resistant cell subpopulations.This relies on a binary vision of structuring phenotypes, exchanges between compartments representing switching between discrete subpopulations, sensitive vs. resistant, as in the SIR models for infectious diseases introduced by Kermack and McKendrick in 1927 [28].The vision presented in this paper is closer to the type of transport equation introduced by McKendrick (alone) in 1926 [39], to which nonlocal terms may be added, with or without diffusion and advection terms.Indeed, there is no reason why functional structuring traits in the population, such as fecundity, viability, motility and plasticity should be discrete.This is the reason why heterogeneity should a priori be represented by continuous traits.It depends on the underlying questions at stake, and on the ways of appreciating such heterogeneity.If the available measures and if the biological questions may be summed up by discrete characters, why not design compartmental ODE models to take account of heterogeneity?However, in general, models structured in continuous traits, yielding integro-differential equations, or partial differential equations with diffusion and advection terms, should be the rule, as heterogeneity in biology is only coarsely represented by discrete variables.

C. Plasticity and tumour microenvironment
As mentioned earlier, plasticity understood as nongenetic instability may be represented by a Laplacian.This is the case in [3], [7], [8], where nongenetic instability or plasticity play a major role in the evolution of the solutions.The influence of evolutionary pressure (or stress) due to the tumour micro-environment is adequately represented in both cases by an advection (or drift) term, the specific effect of which is studied in [9], and interpreted as displacing the expected concentration point of the one-dimensional structure variable of the nonlocal Lotka-Volterra reactiondiffusion under consideration.In [3], plasticity is endowed with the nature of an intrinsic trait tuning the Laplacian and evolving with time as the other (competing) traits.It is noteworthy that in order to represent competing phenotypes, such as fecundity and viability, a dimensionality of the space of phenotypes should be greater than one.If we add Cartesian space to phenotype space (as performed in [36], however with only one dimension in each case), this will make the model even more complex to study.Nevertheless, this is the framework of development of actual solid tumours, in which local (in space) fixation of phenotypes may lead with genetic instability to the emergence of spatially isolated mutations according to a Darwinian-like phylogenetic tree as has been shown in [21].In other words, mixed space and phenotype models should consider firstly concentration of phenotypes, reversible and costly in terms of energy, likely using epigenetic enzymes to face micro-environmental stress [53], in spatially isolated regions (spatial basins of attraction for such concentration), and secondly, under the pressure of microenvironmental stress in these regions, possible emergence of mutations fixing -on a less costly mode than through the action of epigenetic enzymes -a favourable phenotype.

D. Interacting structured cell populations
Interactions of tumour cell populations with surrounding healthy ones may be mutualistic [44] or competitive.Structured models of tumour immune-interactions, for the latter, are in development [27], since the beginning of modern immunotherapies, and they offer many opportunities, including for the development of control methods, to mathematicians.

VI. CONCLUSION AND PERSPECTIVES
In this short review of some mathematical models of tumour growth and its control, focused on structured cell population models, principles of physiologically based modelling for cancer growth, to take into account cell population heterogeneity and plasticity, with examples of models with and without control, have been presented.Clearly, the categories of "cancer" and "healthy" cell populations we have considered in our model with optimal control are very general.The illustrated solution to the above mentioned optimal control problem OCP proposes a qualitative strategy of control that should be adapted to defined clinical settings, in order to be actually useful in practical oncology.Nevertheless, the proposal is on the table, awaiting for therapists to consider it.In the immediate future, modelling tumour-immune interactions with representation of modern immunotherapies involving tumour control by immune checkpoint inhibitors (ICIs, that boost immune cells) is another goal which is presently pursued [27], with the same principle of structuring cell populations according to a continuous phenotype representing their heterogeneity.However, even immunotherapies act by tumour cell kill induced by drugs, and their adverse events, although rare, may be severe, even fatal, as in the case of cardiac or diaphragmatic myopathies.The grail of anticancer treatments might be, in particular by using redifferentiation therapies, to avoid cancer cell kill and its adverse toxic events on healthy tissues.Unfortunately, this has been achieved so far only in very special cases, such as acute promyelocytic leukaemia (APL) or chronic myelogenous leukaemia (CML) with drugs All-TransRetinoic Acid (ATRA) and Imatinib Mesylate, respectively.Progression in the understanding of mechanisms of tissue cohesion in multicellular organisms and their recovery after local impairment in cancer might offer new opportunities for therapeutics, although most of these mechanisms are still widely unknown.