Luis Almeida
email: almeida@ljll.math.upmc.fr

Rebecca H Chisholm
email: chisholm.r@unimelb.edu.au

Jean Clairambault
email: jean.clairambault@inria.fr

Tommaso Lorenzi

Alexander Lorz ¶
email: alexander.lorz.1@kaust.edu.sa

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

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

WHY IS EVOLUTION IMPORTANT IN CANCER AND WHAT MATHEMATICS SHOULD BE USED TO TREAT CANCER? FOCUS ON DRUG RESISTANCE

The clinical question of drug resistance in cancer, our initial motivation to study continuous models of adaptive cell population dynamics, leads naturally and more generally to consider the cancer disease itself from an evolutionary biology viewpoint, a consideration without which even the best targeted therapies will likely most often eventually fail. Among the challenging questions to mathematicians who tackle the task of understanding this disease and optimising its treatment are the representation of phenotypic heterogeneity of cancer cell populations and of their plasticity in response to anticancer drug insults. Such representation can be obtained using phenotype-structured models of healthy and cancer cell populations, and optimal control methods to optimise drug effects, with the perspective to implement them in the therapeutics of cancer, aiming at both avoiding the emergence of drug resistance in tumours and taking into account a constraint of limiting unwanted adverse effects to healthy tissues.

Introduction to mathematical modelling in cancer

Mathematical models of cancer growth and therapy have already known numerous developments and publications in the past twenty years or so. They belong to two general classes: agent-based models, ruled by stochastic rules of growth (for division, death, motion, interactions with the environment) in which the individual agents are cancer cells, and continuous models that rely on ordinary or partial differential equations, sometimes delay differential equations, whose solutions are densities of cancer cell populations. The benefits and limitations of these two respective classes of models, with examples, are discussed, e.g.., in [Chisholm et al. 2016a] 6 . As regards anticancer treatments, the continuous version allows to take advantage of mathematical optimisation and optimal control algorithms that have been designed in this framework, originally in engineering settings. A short review of models designed with this therapeutic control vision is presented for instance in [Billy and Clairambault 2013] 2 . It is sometimes possible to obtain a continuous model starting from an agent-based one by averaging methods; alternatively, one can also develop in parallel the two types of models applied to the same biological problem and compare the predicted behaviour of the modelled cell populations, as e.g., in [Chisholm et al. 2015] 9 , or in a general setting, in [Byrne and Drasdo 2009] 4 .

The goal of such models of cancer growth may be to merely understand the biological phenomenon of cancer growth, by designing accurate models that are all the more relevant to describe a biological reality as they are identified and validated on biological measurements in vitro in culture dishes, in vivo on laboratory animals, or from observations (e.g. radiological images) on humans, to be confronted to theoretical growth curves depending on a priori unknown parameters (the physicist's viewpoint). But it may also be of a different nature, to represent the effects of treatments on tumours, with the aim to optimise them. In the latter case, these effects may be described either by their molecular effects on known drug targets (keeping in mind that precision targeting is often alluring, since drugs may have unpredictable effects on non recognised targets) or by their functional effects on the possible fates of cell populations, namely proliferation, extinction, differentiation or senescence. The respective advantages of these two points of view are also discussed, with examples, in [Chisholm et al. 2016a] 6 . Whatever the chosen point of view, molecular or functional, the goal of these models is here clearly established as understanding and improving the efficacy of anticancer treatments (the physician's viewpoint).

Drug resistance in cancer

The two main pitfalls of cancer therapeutics

Unwanted toxic side effects on healthy cell populations and emergence of resistance to treatments in cancer cell populations are the two main pitfalls of cancer therapeutics in the clinic. Toxicity is always a concern for the clinician, as it limits the tolerable doses of drugs delivered to the patient, who otherwise might see his tumours eradicated, but at the expense of deadly insults to essential organs or functions (haematopoiesis, digestion, skin covering, liver function, heart function). It has lately been proposed that instead of delivering for short periods of time the maximum tolerated dose (MTD), it might be as efficient to deliver small drug doses (the socalled metronomic strategy), thus minimising toxicity, with as good results on the cancer cell population. To what extent is the immune system involved in the efficacy of this new way of designing delivery schedules is not completely clear and might depend on the anticancer drug in use [Zitvogel et al. 2008, 2011] 39,40 , however metronomic therapies certainly challenge the MTD strategy in both limitation of toxicity and improvement of efficacy. Note that an initial interpretation of the success of metronomic therapies was more mechanic, postulating that too high amounts of drugs destroy the blood vessels that bring the drugs to the tumour [Benzekry and Hahnfeldt 2013, Pasquier et al. 2010] 1,28 . It is not exclusive of the immunogenic explanation [Zitvogel et al. 2008, 2011] 39,40 , which proposes that giving small drug doses may reveal hidden (internalised) cancer antigens by shattering a small number of cancer cells, enough to trigger an efficient immune response towards the whole cancer cell population. Both explanations still remain to be more biologically documented -especially with respect to the immune response -, mathematically modelled and tried in clinical settings; however, they address the question of toxicity in an apparently paradoxical way ("more is not necessarily better") that is a challenge for modellers.

Drug resistance, the other major pitfall of cancer therapeutics, is a treatment efficacy limitation of another nature; it may be defined as adaptation of the target cancer cell populations to the hostile environment created by the drug. Resistance to treatments in cancer cell populations, insofar as it is not constitutive of organisms therapies apply to, but secondary, i.e., induced by treatments as a stress response. In many cases (in fact, in most cases), treatments that show remarkable initial efficacy by drastically shrinking tumours see their response decrease with time, until they become totally inefficient as tumours regrow. Furthermore, the newly growing can-cer cell populations, that have become resistant to the drug in use, are out of reach for this therapy, and often for others that have not been employed (multi-drug resistance). At the molecular scale, different mechanisms have been identified, such as overexpression of drug efflux pumps (ABC transporters [Gottesman et al. 2002] 16 , such as the P-glycoprotein, also known as MDR1, or ABCB1), of intracellular drug processing enzymes or of DNA repair enzymes, and it has been proposed to combine cytotoxic drugs with inhibitors of these mechanisms, unfortunately eventually to no avail. As mentioned above, the molecular point of view in pharmacological treatments in principle offers a satisfying framework to perform cancer treatment optimisation, but so-called targeted therapies (i.e., that target intracellular molecular pathways), with a few exceptions, result in disappointing outcomes (see, e.g., [Gillies et al. 2012, Goldman et al. 2017a] 13,14 ).

From the single cancer cell to cancer cell populations

Indeed, these treatments share the same flaw, which is that they focus on a given molecular target (or on several molecular targets), considering cancer as the disease of the same single cell extended to large quantities, instead of taking into account the population of cancer cells in its diversity, which might offer a key explanation of their failure [Goldman et al. 2017b] 15 . Such population diversity (or heterogeneity) is not necessarily of genetic nature but linked to epigenetic changes in the chromatin, thus reversible [Sharma et al. 2010] 34 , at least on the initiation of drug resistance (mutations can come later to irreversibly establish resistance in a subclone of the cell population), and may result in differently expressed phenotypes in different cells, potentially inducing different resistance mechanisms as responses to cytotoxic stress in a population of cells bearing all the same genotype.

Introducing the population of cancer cells (indeed, the actual target of anticancer treatments) naturally sets the scenery for Darwinian evolution of cells exposed to anticancer drugs seen as an environmental selection pressure, as will be developed in the next section. This viewpoint, introduced in theoretical ecology for quite a long time already, is rather new in biology and medicine (where it has given rise to the new field of Darwinian medicine), however, does not allow to decide whether the selection is of pure Darwinian nature (selection of the fittest, cells that were already present in the population before exposure to the drug) or may involve a part of Lamarckian adaptation (no resistant cells initially present, but stochastic triggering of resistance mechanisms in a few cells for which the response to stress happens to be well adapted to resist the cytotoxic effects of the drug in use). This alternative, discussed in a mathematical setting in [Chisholm et al. 2015] 9 was already the object of the biological experiment by Luria and Delbrück [Luria and Delbrück 1943] 27 , concluding to sheer Darwinian selection. However, Luria and Delbrück's experiment was performed not on cancer cells exposed to drugs, but on bacteria exposed to bacterium-eating viruses (phages), while human and animal cancer cells bear a genome -and epigenome, i.e., chromatin (histones) -than is by far richer than the bacterial genome, which in our case does not allow to conclude. Nevertheless, the cell population point of view clearly opens new ways to understand and overcome drug resistance in cancer.

Cancer as evolutionary disease

Evolution of multicellularity and cancer

Darwinian evolution (together with possible Lamarckian adaptation) of cancer -and healthy -cell populations (but healthy cell populations are in principle well controlled as regards their possibilities of phenotype evolution) must of course be considered on the short-time level of a human life or disease, but the much larger time of evolution in the course of billions of years, from unicellular organisms towards the organised and coherent forms of multicellularity represented by present animals and plants, may also shed light on cancer as evolutionary disease. Cancer is a disease of multicellular organisms, that may be defined as loss of coherence between tissues due to loss of coherence control by those genes that have been essential in the evolution towards multicellularity. In [Davies and Lineweaver 2011] 10 , it is advocated that the genes that are altered in cancer are precisely the ones that have been employed by evolution to design multicellular organisms. Indeed, evolution proceeds, as stated by Nobel prize laureate François Jacob in [Jacob 1977] 20 , by tinkering, i.e., it proceeds by trials and errors taking advantage of any existing material, and, as regards multicellularity, such tinkering may result in localised (in organs and physiological functions elicited by corresponding genes) fragilities, that secondarily, under environmental pressure, may be caught off guard and result in localised cancers. Such loss of coherence control, unmasking a pluripotent phenotype that is also named plasticity, may in particular be seen in the process of de-differentiation of cancer cells, i.e., adoption of a pluripotential phenotype (eventually yielding the so-called 'cancer stem cell', whose existence is likely to be transient [Li and Laterra 2012] [START_REF] Li | [END_REF] ) making the cells that bear it, as endowed with a rich panel of non-repressed genes, able to develop a wide variety of responses to cytotoxic stress. The involvement of such failed multicellularity (i.e., unpreserved normal differentiation) control genes in revealing an ancient 'toolkit' of prexisting adaptations [Davies and Lineweaver 2011] 10 still remains to be documented, but it certainly offers new ways of considering cancer as an evolutionary disease and drug resistance in cancer as an evolutionary -and adaptive -mechanism.

Heterogeneity and plasticity in cancer cell populations

Heterogeneity in cancer cell populations has been documented in advanced solid tumours as of genetic nature, with evidence of multiple branched mutations [Gerlinger et al. 2012] 12 , but, as mentioned above, it may also consist of sheer epigenetic and reversible modifications [Sharma et al. 2010] 34 linked to enzymatic activities located on the chromatin, i.e., without mutations in the genome. However, as recently shown in [Shaffer et al. 2017] 33 , such fast epigenetic, non genetic, reprogramming of a sparse subpopulation of cancer cells may eventually result in a stably resistant state.

Another look at heterogeneity induced in cancer cell populations by exposure to cytotoxic drugs is presented in [Wu et al. 2015] [START_REF] Wu | Proc. Nat. Acad. Sci. USA[END_REF] . In this article, it is proposed that so-called cold genes, that have been identified as expressed in the genome of cancer cells (multiple myeloma cells) have a very ancient origin, being conserved without changes throughout evolution from unicellularity, and may be responsible for stress response in extremely hostile and unpredictable situations (resulting from events comparable, mutatis mutandis, to the impact on animal life of the meteorite that 60 million years ago fell on Earth -creating the Chicxulub crater in Yucatan -, subsequently putting an end to the dominance of dinosaurs), by possibly launching secondary expression of various resistance mechanisms. In this respect, these very ancient 'cold genes', elaborated in a remote past of our planet, when conditions of life were different from the present (UV radiation, acidity, low oxygen concentration in the oceans and in the atmosphere), might be the genetic toolkit of preexisting adaptations mentioned above, or part of it.

The variety of resistance mechanisms developed by cancer cell populations exposed to letal doses of cytotoxic drugs has been related to what is called bet hedging in theoretical ecology. The term 'bet hedging' is used to qualify behaviour relying on an ensemble of traits that make a population of living individuals adaptable to an unpredictable environment, using a so-called 'risk-spreading strategy' , that at the scale of the population, may result in keeping safe only a small part of it, but a part that will be able to reconstitute the whole population, with preservation of its common genome, after such adaptation to the new environment [Philippi and Seger 1989] 29 . Bet hedging in tumours has been proposed as a stochastic 'cancer strategy' in [Brutovsky and Horvath 2013] 3 . It is also presented as "an ultimate explanation of intra-tumour heterogeneity" in chapter XVII of the book [Ujvari, Roche and Thomas 2017] [START_REF]Ecology and evolution of cancer[END_REF] .

Plasticity, mentioned above about de-differentiation of cancer cells and the transient state of cancer stem cell, may be evidenced at the level of the single cell (derepression of genes that must be epigenetically repressed in physiology to produce the differentiation that yields about 200 different functional cell types in the human organism), but also at the level of the cell population, since the spreading of such pluripotent cells makes the population adaptable to environmental changes (plastic), possibly by using expression of cold genes in a tiny subpopulation and stochastic (or distributed) bet hedging of resistance phenotypes.

The plasticity -physiologically normal in highly undifferentiated cell states, close to stem cells, but totally pathological in cell populations for which a defined terminal physiological function exists -of the epigenetic landscape of a given human genome, as metaphorically proposed in [Waddington 1957] [START_REF] Waddington | The strategy of genes[END_REF] , recently revised from a systems biology viewpoint by Sui Huang, see, e.g., [Huang 2011, 2013] 17,18 , provides another approach to plasticity and evolution of cancer, that has been for instance exploited to study lineage commitment in haematopoiesis by using bifurcation analysis of an ordinary differential equation model [Huang et al. 2007] 19 .

Continuous mathematical models

Phenotype-structured mathematical models

The modelling framework of adaptive dynamics we present here is more likely to correspond biologically to epigenetic modifications rather than to genetic mutations, as the evolution in phenotype is in this mathematical setting always reversible (not to mention that eventual induction of emergent resistant cell clones due to mutations under drug pressure is never to be excluded in the long run). From the biologist's point of view, we study phenotypically heterogeneous, but genetically homogeneous, cancer cell populations under stress by drugs.

The models considered here are all based on the so-called logistic ODE model, which we recall here. It is given by the equation

dN dt = (r -dN ) N,
which describes the time-evolution of the number of individuals N (t), starting from a prescribed initial condition N 0 . Coefficient r denotes the net selection rate of the individuals, namely the difference between their proliferation and death rates, while the logistic term dN stands for an added death rate proportional to the number of individuals.

The underlying assumption is that competition for nutrients and space inside the population does not allow for exponential unconstrained proliferation. Mathematically, it is indeed true that if N 0 ≤ r d , then N (t) converges increasingly toward the carrying capacity r d . Let us now introduce a basic phenotype-structured model, where the quantity of interest is a number of cells n(t, x) at time t > 0, and of phenotype x ∈ (0, 1) standing for the resistance to a given drug. We stress that this phenotype is taken to be continuous, because, as already mentioned, it can be correlated to biological characteristics which themselves are continuous. Here, (0, 1) is taken for simplicity but multi-dimensional phenotypes can of course be considered.

The model reads

∂n ∂t (t, x) = (r(x) -d(x)ρ(t)) n(t, x),
where ρ(t) := 1 0 n(t, x) dx is the total number of individuals at time t, starting from some initial condition n 0 (•). As before, r(x) is the net proliferation rate of cells of cells of phenotype x, while d(x)ρ(t) is the natural extension of the previous logistic term. Note that more general logistic death terms through a Kernel K can be considered, in the form 1 0 K(x, y)n(t, y) dy. Model (4.1) is characterised by two main phenomena: convergence and concentration. The first one means convergence of ρ(t) towards max r d , and concentration of the density n on the set of phenotypes where r d reaches its maximum, namely arg max r d . This is why this class of models is extensively used in adaptive dynamics to model selection: only some fitter phenotypes remain, which corresponds mathematically to the convergence of n(t, •) to a sum of Dirac masses located on the set arg max r d . This modelling framework also extends to several populations, in which case the competition between the populations are modelled through Lotka-Volterra-like terms. Let us introduce a model of two interacting popula-tions, which will be further developed in the next section with the modelling of chemotherapy as control terms. It is concerned with two densities of healthy and cancer cells n H (t, x) and n C (t, x) respectively, where x is again a continuous phenotype describing the resistance to a given drug. The model is given by

∂n H ∂t = r H (x) -d H (x) (a HH ρ H (t) + a HC ρ C (t)) n H (t, x), (1)
∂n C ∂t = r C (x) -d C (x) (a CC ρ C (t) + a CH ρ H (t)) n C (t, x), (2)
where, as before, ρ 31 that the behaviour of ( 1) is again convergence and concentration, where the asymptotic values of ρ H , ρ C and the sets on which n H , n C concentrate can also be explicitly computed.

H (t) = 1 0 n H (t, x) dx, ρ C (t) = 1 0 n C (t,

Optimal control for anticancer therapeutics

Optimal control methods (reviewed in [Trélat 2005] [START_REF] Trélat | Contrôle optimal, Vuibert[END_REF] ) applied to models of cancer therapeutics using systems of ordinary differential equations [Carrère 2017; Ledżewicz and Schättler 2006, 2014] 5,21,22 or of partial differential equations [Pouchol et al. 2016] 31 are the appropriate tool to theoretically optimise cancer therapeutics, in particular by taking into account the inevitable emergence of drug resistance in cancer cell populations.

The built-in targets for theoretical therapeutic control that are present in the phenotype-structured PDE models we advocate here are not supposed to represent well-defined molecular effects of the drugs in use, but rather functional effects, i.e., related to cell death (cytotoxic drugs), or to proliferation in the sense of slowing down the cell division cycle without killing cells (cytostatic drugs). We propose that cell life-threatening drugs (cytotoxics) induce by far more resistance in the highly plastic cancer cell populations than drugs that only limit their growth (cytostatics), and that a rational combination of the two classes of drugs -and possibly others, adding relevant targets to the model -may be optimised to propose therapeutic control strategies to avoid the emergence of drug resistance in tumours.

We address this optimal control problem in the context of two populations, healthy and cancer as in the model is given by ( 1), now complemented with two types of drugs of infusion rates u 1 and u 2 for cytotoxic and cytostatic drugs, respectively. The resistance phenotype x they are endowed with is defined with respect to the cytotoxic drug pressure, and is taken to range from sensitiveness (x = 0) to resistance (x = 1). The controlled model thus reads

       ∂n H ∂t (t, x) = r H (x) 1+α H u2(t) -d H (x)I H (t) -µ C (x)u 1 (t) n H (t, x), ∂n C ∂t (t, x) = r C (x) 1+α C u2(t) -d C (x)I C (t) -µ H (x)u 1 (t) n C (t, x), (3)
On a fixed therapeutic time-window [0, T ], the optimal control problem is to choose the controls u 1 and u 2 so as to minimise the number of cancer cells ρ C (T ), while satisfying the three following constraints.

• remaining under maximum tolerated doses:

0 ≤ u 1 (t) ≤ u max 1 , 0 ≤ u 2 (t) ≤ u max 2 ,
• avoiding the emergence of too big a tumour :

ρ H (t) ρ H (t)+ρ C (t) ≥ θ HC ,
• limiting unwanted adverse effects to the healthy cell population:

ρ H (t) ≥ θ H ρ 0 H .
This optimal control problem is motivated by the inefficacy of using constant high doses of drugs, a strategy which on the long run violates the last two constraints. This is indeed what is observed in the simulation presented in Figure 1: although the tumour size first decreases, it is at the expense of the cancer cell density concentrating on a resistant phenotype. The treatment becomes inefficient and relapse occurs.

In [Pouchol et al. 2016] 31 , the previously defined optimal control problem is analysed both numerically and theoretically. As the time T increases, it is found that the optimal control strategy becomes increasingly close to a two-phase strategy.

• The first phase is long, and only constant low doses are given, so as to saturate the second constraint. At the end of this first phase, the drug pressure has been low enough to ensure that the cancer cell density has concentrated on a sensitive phenotype. • The second phase is short and starts with maximum tolerated doses for both drugs, leading to a quick decrease of both cell numbers be- cause they are efficient on a sensitive cancer cell population. Once the third constraint (on the heathy cell density) has been reached, cytostatic drugs switch to some intermediate value (which can be computed in feedback form) which allows for a further decrease of the tumour size while keeping the healthy cell number at its lower bound.

A numerical simulation of the optimal strategy is presented in Figure 2 below.

For a pratical implementation of the previous strategy, it is natural to repeat it in a quasi-periodic manner. One can hope that after enough cycles, the tumour will be eradicated, or at least made chronic. In order to decide when to switch from the second short phase to another cycle with a long first phase, one must identify markers for resistance. Indeed, as long as constant low doses do not violate the second constraint on the relative tumour size, they must be given to ensure that the (assumed to be plastic) tumour is becoming sensitive to the treatment again. The switch to the second phase can be led as soon as the markers indicate that the tumour has become sensitive enough again. Finally, if the healthy cells tissue is too damaged (namely the third constraint saturates), one can hope to still let the tumour decrease with a properly chosen cytotoxic drug infusion. When this is no longer possible, one must switch the the long first phase.

5. Future tracks in modelling for cancer therapeutics

Beyond present models to optimise cancer therapeutics

The models of adaptive population dynamics that we have presented here, with their built-in targets for control, rely on a nonlocal Lotka-Volterra vision of cell-cell population competition. This point of view could be extended to other modes of interaction, which could be mutualistic or predator-prey like, and to an arbitrary number of cell populations. For an analysis of a mutualistic integro-differential 2x2 system modelling interactions between breast cancer cells and their supporting stroma (adipocytes), we refer to [Pouchol 2015] 30 , whereas a more general Lotka-Volterra-like model for N populations is analysed in [Pouchol and Trélat] 32 . The inferred asymptotic behaviour is again convergence and concentration.

Let us now come to extensions of the integro-differential setting by considering the basis model (4.1). Recall first that more general logistic interaction terms can be considered. The other natural extension is to model epimutations (occurring on the relevant time-scale, which is here that of a tumour). They can either be modelled by a Laplacian, leading to

∂n ∂t (t, x) = (r(x) -d(x)ρ(t)) n(t, x) + β∆n(t, x),
with Neumann boundary conditions, or more generally through a mutation Kernel. Note that both modelling are linked at the limit through a proper rescaling of the Kernel, as explained in [Perthame 2007] ? . A complementary advection term can be added, accounting for cells actively adapting to their environment, seeking for phenotype changes that make them fitter. These can be seen as stress-induced epimutations and with them the model becomes

∂n ∂t (t, x) + ∂ ∂x (v(x)n(t, x)) = (r(x) -d(x)ρ(t)) n(t, x) + β∆n(t, x).
Note that in [Chisholm et al.] 9 , the advection happens to be compulsory to observe quick enough dynamics to fit those obtained in the experiments presented in [Sharma et al.] 34 .

A further advantage of these more general PDE models is that they are able to represent possible asymptotic coexistence of phenotypes, which is not the typical output of the integro-differential models.

A final possible extension worth-mentioning is the addition of a space variable, since it is believed that the heterogeneity of a tumour varies from its periphery to its center. This is also particularly relevant in view of optimal control through chemotherapy since drugs will efficientlly access the outer rim of the tumour but less its core. For possible cancer models taking both phenotype and space into account, we refer to [Lorz et al. 2015] ? .

Need for models with a larger evolutionary perspective

From the biological part of this article, it clearly appears that the abovementioned models, sophisticated though they may be, are not enough to study in mathematical settings the evolution of multicellularity and its intrinsic failure, namely cancer, nor are they presently enough to design optimised therapeutic strategies that can overcome drug resistance in cancer. Open questions to biologists remain if one wants to make available a framework within which mathematical modelling may be designed. What are the genes that must be silenced in physiology and are re-expressed in de-differentiated cancer cells? What are the observable links between genes that are known to be essential for multicellularity and what are the genes that are altered in cancer (the same, following [Davies and Lineweaver 2011] 10 ])? What models to study physiological coherence between tissues in the same organism, i.e., what sticks together in harmony the 200 different cell types of a human organism? What part of the genome bears the so-called cold genes, what part the individual signature of an organism that is transmitted throughout differentiation (the 'self', certainly to be related with the major histocompatibility complex, MHC), and what part the genes that are normally sequentially silenced in the history of differentiations? As regards mechanisms of drug resistance, what part of launching in a cancer cell population is deterministic (triggering cold genes) and what part is stochastic? And such list of open questions is not intended to be comprehensive.

Conclusion: a challenging new field for mathematicians

In this short description of cancer as evolutionary disease, focusing on the question of drug resistance and its possible overcoming by optimised strategies in the clinic, we have presented what has been recently developed in the framework of mathematical modelling, that is, adaptive dynamics of cell populations represented by phenotype-structured models relying on partial differential equations, together with optimal control methods to guide their asymptotic behaviour. We have also proposed immediate tracks for future extensions of these existing models, and only sketched the scenery for future mathematical models that still lack biological answers to guide their design. We are nevertheless confident in the fast progress of cancer biology to help mathematicians design models that can be helpful in prevention, prediction and control of cancer in the clinic, provided that the right questions are posed, mathematically challenged and experimentally tackled.

Figure 1 .

 1 Figure 1. Simulation with u 1 (t) = Cst = 3.5 and u 2 (t) = Cst = 2, in time T = 10. Here ρ CS (t) = 1 0 (1 -x)n C (t, x) dx is a measure of the number of sensitive cells in the cancer cell population.

Figure 2 .

 2 Figure 2. Simulation of the optimal control problem for T = 60.

  x) dx. The logistic terms now incorporate an intraspecific competition term weighted by coefficients a HC and a CH . Because cells belong to different ecological niches, it is quite natural tu assume a HC < a HH , a CH < a CC .Under hypothesis (4.1), it is proved in[Pouchol et al. 2016]


Simone Baldi

Ridho Muhammad
email: muhammadridhorosa@tudelft.nl

Paolo Rosa

Frasca

Muhammad Ridho

Paolo Frasca
email: paolo.frasca@gipsa-lab.fr

Adaptive state-feedback synchronization with distributed input: the cyclic case

Keywords: Adaptive synchronization, heterogeneous uncertain agents, adaptive control

    Adaptive state-feedback synchronization with distributed input: the cyclic case

INTRODUCTION

In recent years, cooperative control of multiagent systems has received increasing attention, due to its impact in formation flying, smart energy, traffic, and other areas [START_REF] Ren | Information consensus in multivehicle cooperative control[END_REF]; [START_REF] Bullo | Gossip coverage control for robotic networks: Dynamical systems on the space of partitions[END_REF]). An important problem in cooperative control is to achieve in a distributed way (i.e. using local information) a common behavior for the entire network: this is the so-called synchronization problem [START_REF] Dorfler | Synchronization in complex networks of phase oscillators: A survey[END_REF]; [START_REF] Turci | Adaptive pinning control: A review of the fully decentralized strategy and its extensions[END_REF]; [START_REF] Gibson | Adaptation and synchronization over a network: Asymptotic error convergence and pinning[END_REF]; [START_REF] Casadei | Multi-pattern output consensus in networks of heterogeneous nonlinear agents with uncertain leader: a nonlinear regression approach[END_REF]). Synchronization has been studied for: uncertain but homogeneous agents [START_REF] Li | Distributed adaptive consensus and output tracking of unknown linear systems on directed graphs[END_REF]; [START_REF] Ding | Distributed adaptive consensus control of nonlinear output-feedback systems on directed graphs[END_REF]), or heterogeneous agents with limited uncertainty [START_REF] Seyboth | Cooperative control of linear multi-agent systems via distributed output regulation and transient synchronization[END_REF][START_REF] Seyboth | On robust synchronization of heterogeneous linear multi-agent systems with static couplings[END_REF]; [START_REF] Li | Distributed robust consensus control of multi-agent systems with heterogeneous matching uncertainties[END_REF]; [START_REF] Mei | Distributed consensus of second-order multi-agent systems with heterogeneous unknown inertias and control gains under a directed graph[END_REF]). It results that synchronization for agents that are concurrently heterogeneous and uncertain is still a major problem. Recently, to handle heterogeneity and uncertainty, it has been proposed to formulate the synchronization problem as a special model reference adaptive control (MRAC) in which each agents tries to converge to the model defined by its neighbors (Baldi and Frasca (2018); Harfouch et al. (2017a)). This formulation is based on 'feedback matching gains' (used to match each agent to the reference model, or leader) and 'coupling matching gains' (used to match each agent with its neighbors). Adaptive laws for both feedback and coupling gains are derived via Lyapunov analysis. The distributed MRAC exploits a communication setting in which the input is communicated among neighbors (distributed input approach). This is alternative to another popular approach to synchronization, the distributed observer approach [START_REF] Cai | The adaptive distributed observer approach to the cooperative output regulation of linear multi-agent systems[END_REF]; [START_REF] Lu | Cooperative output regulation of linear multi-agent systems by a novel distributed dynamic compensator[END_REF]; [START_REF] Baldi | Cooperative output regulation of heterogeneous unknown systems via passification-based adaptation[END_REF]) where, in place of the input, an observation of the leader state is communicated among neighbors. Both the distributed input and the distributed observer approaches include a feedforward action and need to communicate auxiliary variables to the neighbors (inputs and observations, respectively) to reconstruct the reference signal. Despite these similarities, the distributed observer scheme can handle cyclic graphs and undirected graphs, while the distributed input with MRAC can be applied only to acyclic directed graphs. For this reason it finds main application in platooning, where no cyclic communication occurs [START_REF] Harfouch | An adaptive switched control approach to heterogeneous platooning with inter-vehicle communication losses[END_REF]. Even if the distributed observer approach can be used in a larger number of cases, the distributed input approach is relevant because the dimension of the input vector is typically smaller than the dimension of the leader state vector, and thus communication with distributed inputs is less cumbersome. While hierarchical architectures have been proposed to remove cycles in the distributed input approach [START_REF] Wang | Hierarchical decomposition based consensus tracking for uncertain interconnected systems via distributed adaptive output feedback control[END_REF]), the open question that motivates this work is: is it possible to design MRAC algorithms based on distributed input with the ability to handle cyclic and undirected graphs? This work gives a positive answer by showing that the same MRAC parametrization derived for the acyclic graph case can be extended to cyclic and undirected graphs. Notation: The transpose of a matrix or of a vector is indicated with X ′ and x ′ respectively. A directed graph (digraph) is indicated with the pair (N , E), where N is a finite nonempty set of nodes, and E ∈ N × N is a set of ordered pair of nodes, called edges. The adjacency matrix A = [a ij ] of an unweighted digraph is defined as a ii = 0 and a ij > 0 if (j, i) ∈ E, where i = j.

THE ACYCLIC CASE

Fig. 1. A simple acyclic network

To recall the synchronization results in the acyclic case, let us consider the network in Fig. 1. Three agents, denoted with indices 1, 2 and 3, have uncertain dynamics

ẋ1 = A 1 x 1 + b 1 u 1 ẋ2 = A 2 x 2 + b 2 u 2 ẋ3 = A 3 x 3 + b 3 u 3 (1)
where The synchronization task is achieved when x 1 , x 2 , x 3 → x m for t → ∞. Being the system matrices in (1) unknown, the synchronization task has to be achieved in an adaptive fashion. In order to have a well-posed adaptive control problem, the following assumptions should be verified. Assumption 1. [Feedback matching conditions] There ex-

x 1 , x 2 , x 3 ∈ R n is the state, u 1 , u 2 , u 3 ∈ R
ist vectors k * 1 , k * 2 , k * 3 and scalars l * 1 , l * 2 , l * 3 such that A m = A 1 + b 1 k * ′ 1 , b m = b 1 l * 1 A m = A 2 + b 2 k * ′ 2 , b m = b 2 l * 2 A m = A 3 + b 3 k * ′ 3 , b m = b 3 l * 3 . (3)
Assumption 2. The signs of the input vector fields, i.e. the signs of l * 1 , l * 2 , l * 3 , are known. Assumptions 1 and 2 are classical conditions mutuated from MRAC [START_REF] Tao | Adaptive Control Design and Analysis[END_REF]; [START_REF] Ioannou | Robust Adaptive Control[END_REF]). We deal with the single-input case, although extension to the multi-input is possible following the related multivariable adaptive control theory. A consequence of Assumption 1 is the existence of coupling gains between neighboring agents.

A 1 = A 2 + b 2 k * ′ 21 , b 1 = b 2 l * 21 A 1 = A 3 + b 3 k * ′ 31 , b 1 = b 3 l * 31 A 2 = A 3 + b 3 k * ′ 32 , b 2 = b 3 l * 32 . (4)
Proof. To derive (4), we find from ( 3)

b 1 = b 2 l * 2 l * 1 , A 1 -A 2 = b 2 k * ′ 2 -b 1 k * ′ 1 (5) which gives k * ′ 21 = k * ′ 2 - l * 2 l * 1 k * ′ 1 and l * 21 = l * 2 l * 1 . Similar calculations hold for k * 31 , k * 32 , l * 31 , l * 32 .
The synchronization of agent 1 to the reference model is the well-known model reference adaptive control (Tao, 2003, Chap. 4): it amounts to the controller

u 1 (t) = k ′ 1 (t)x 1 (t) + l 1 (t)r(t) (6) with k 1 , l 1 the estimates of k * 1 , l * 1
, and to the adaptive laws

k′ 1 = -sgn(l * 1 )γ k b ′ m P e 1 x ′ 1 l1 = -sgn(l * 1 )γ l b ′ m P e 1 r (7
) where e 1 = x 1x m , the scalars γ k , γ l > 0 are adaptive gains, and P is a positive definite matrix satisfying

P A m + A ′ m P = -Q, Q > 0. (8) Proving synchronization exploits the Lyapunov function V 1 (e 1 , k1 , l1 ) = e ′ 1 P e 1 + tr k′ i ki γ k |l * i | + l2 i γ l |l * i | (9)
and the error dynamics ė1 = A m e 21 + b 1 ( k′ 1 x 1 + l1 r).

(10) The details are well known, cf. (Tao, 2003, Chap. 4).

The synchronization of agent 2 to agent 1 should avoid the use of r. This is possible via the controller 11) and the adaptive laws we can calculate the time derivative of ( 14) along ( 13)

u 2 (t) = k ′ 21 (t)x 1 (t)+k ′ 2 (t)(x 2 (t)-x 1 (t))+l 21 (t)u 1 (t) (
k′ 21 = -sgn(l * 2 )γ k b ′ m P (x 2 -x 1 )x ′ 1 k′ 2 = -sgn(l * 2 )γ k b ′ m P (x 2 -x 1 )(x 2 -x 1 ) ′ l21 = -sgn(l * 2 )γ l b ′ m P (x 2 -x 1 )
V21 = -e ′ 21 Qe 21 + 2(sgn(l * 2 )b ′ m P e 21 x ′ 1 + γ -1 k k′ 21 ) k′ 21 |l * 2 | +2(sgn(l * 2 )b ′ m P e 21 e ′ 21 + γ -1 k k′ 2 ) k2 |l * 2 | +2(sgn(l * 2 )b ′ m P e 21 u 1 + γ -1 l l21 ) l′ 21 |l * 2 | = -e ′
21 Qe 21 (15) where we have used (12). Using standard Lyapunov arguments and Barbalat's lemma we can show V21 → 0 as t → ∞ and hence e 21 → 0.

To synchronize agent 3 to agents 1 and 2, let us derive the dynamics of the error e 31 = x 3x 1 and e 32 = x 3 - 16) which leads us to the controller

x 2 ė31 = A m e 31 + b 3 (u 3 -k * ′ 31 x 1 -k * ′ 3 e 31 -l * 31 u 1 ) ė32 = A m e 32 + b 3 (u 3 -k * ′ 32 x 2 -k * ′ 3 e 32 -l * 32 u 2 ) (
u 3 (t) = k ′ 31 (t) x 1 (t) 2 + k ′ 32 (t) x 2 (t) 2 + k ′ 3 (t)
k′ 31 k31 γ k |l * 3 | + tr k′ 32 k32 γ k |l * 3 | +tr k′ 3 k3 γ k |l * 3 | + l2 31 γ l |l * 3 | + l2 32 γ l |l * 3 | . ( 19
)
It is possible to verify V321 = -e ′ 321 Qe 321 and e 321 → 0 using similar Lyapunov arguments as before. This shows that the state of agent 3 converges to the average of the states of agents 1 and 2. Since the states of agents 1 and 2 converge to the state x m of the reference model, then the state of agent 3 will converge to x m as well (cf. [START_REF] Rosa | Synchronization of uncertain heterogeneous agents: an adaptive virtual model reference approach[END_REF] for the full details). Overall synchronization can be proved via the Lyapunov function

V 1 + V 21 + V 321 .
Remark 1. The adaptation laws (12) and (18) reminds the setting of systems stabilizing each other through adaptation [START_REF] Narendra | Unstable systems stabilizing each other through adaptation[END_REF]), with the peculiar difference that the directed path to the leader makes our adaptation always stable. On the other hand, [START_REF] Narendra | Unstable systems stabilizing each other through adaptation -part ii[END_REF]) discusses instability due to the absence of such leader.

General acyclic case

Extending from Fig. 1, let us consider a set of N agents ẋi = A i x i + b i u i , i ∈ {1, . . . , N } (20) where agent 1 is the one that can access the reference r in (2). Assumptions 1 and 2 are generalized to

A m = A i + b i k * ′ i , b m = b i l * i (21
) with known signs of l * i . Remark 2. Similarly to Proposition 1, one can verify the existence, for every pair of agents (i, j), of a constant vector k * ji and a scalar l * ji such that

A i = A j + b j k * ′ ji , b j = b 2 l * ji . ( 22
)
For convenience of notation, let us rewrite (2) as ẋ0

= A 0 + x 0 + b 0 u 0 (23) with x 0 = x m , A m = A 0 + b 0 k * ′ 0 , b 0 l * 10 = b m , u 0 = k * ′ 0 x 0 + l * 0 r,
where k * 0 , l * 0 are known and do not need to be estimated. This gives the controller (equivalent to (6)) Fig. 2. Simple undirected (left) and cyclic (right) networks

u 1 (t) = k ′ 10 (t)x 1 (t) + k ′ 1 (t)(x 1 (t) -x 0 (t)) + l 10 (t)u 1 (t).
(24) Under the following assumption a synchronization result holds. Assumption 3. The directed communication graph is acyclic. In addition, the graph contains a directed spanning tree with the leader as the root node. Theorem 1. Under Assumptions 1-3, consider the linear systems ( 20), with reference model ( 23), controllers

u j (t) = N i=0 a ij k ′ ji (t)x i (t) N i=0 a ij + k ′ j (t) N i=0 a ij (x j (t) -x i (t)) N i=0 a ij + N i=0 a ij l ji (t)u i (t) N i=0 a ij (25)
with the index i = 0 used for the reference model (i.e. a j0 = 0 only for the root node), and update laws

k′ ji = -sgn(l * j )γ k b ′ m P N i=0 a ij e ji x ′ i k′ j = -sgn(l * j )γ k b ′ m P N i=0 a ij e ji N i=0 a ij e ji ′ lji = -sgn(l * j )γ l b ′ m P N i=0 a ij e ji u i (26)
where e ji = x jx i , and k ji , k i , l ji are the estimates of k * ji , k * i , l * ji respectively. Then, all closed-loop signals are bounded and e i = x ix m , e ji = x jx i , with i, j such that a ij = 0, converge asymptotically to zero. Proof 1. The proof uses the Lyapunov function

V = N j=1 N i=0 a ij e ji ′ P N i=0 a ij e ji + N j=1 tr k′ j kj γ k l * j + N j=1 N i=0 a ij tr k′ ji kji γ k l * j + N j=1 N i=0 a ij l2 ji γ l l * j (27)
Stability tools are similar as before and left to the reader.

THE CYCLIC AND UNDIRECTED CASE

To understand the effect of cycles and undirected links, let us consider the undirected network of Fig. 2 (left). Let us assume we can calculate the inputs using the same method of Theorem 1. The ideal control actions are

2u 1 = k * ′ 1 x 1 + l * 1 r + k * ′ 1 (x 1 -x 2 ) + k * ′ 12 x 2 + l * 12 u 2 u 2 = k * ′ 2 (x 2 -x 1 ) + k * ′ 21 x 1 + l * 21 u 1 . ( 28
)
In order to unequivocally determine u 1 and u 2 , the following equation should have a unique solution

2 -l * 12 -l * 21 1 u 1 u 2 = k * ′ 1 (2x 1 -x 2 ) + l * 1 r + k * ′ 12 x 2 k * ′ 2 (x 2 -x 1 ) + k * ′ 21 x 1 .
From Proposition 1 we have l * 12 l * 21 = 1, so the determinant of the square matrix above is 2l * 12 l * 21 = 1, and the ideal inputs u 1 and u 2 are well defined. In addition, the synchronization error dynamics with the ideal gains is

ė1 + ė12 = 2A 1 x 1 + b 1 k * ′ 1 x 1 + b 1 l * 1 r -A m x m -b m r -b 1 k * ′ 1 x 2 + k * ′ 12 x 2 + l * 12 u 2 -A 2 x 2 -b 2 u 2 ė21 = A 2 x 2 + b 2 k * ′ 2 (x 2 -x 1 ) + b 2 k * ′ 21 x 1 + b 2 l * 21 u 1 -A 1 x 1 -b 1 u 1
(29) which leads to ė1 + ė12 = A m (e 1 + e 12 ), ė21 = A m e 21 .

(30)

Let us now consider the directed cyclic network of Fig. 2 (right) and calculate the inputs using the method of Theorem 1

2u 1 = k * ′ 1 x 1 + l * 1 r + k * ′ 1 (x 1 -x 3 ) + k * ′ 13 x 3 + l * 13 u 3 u 2 = k * ′ 2 (x 2 -x 1 ) + k * ′ 21 x 1 + l * 21 u 1 u 3 = k * ′ 3 (x 3 -x 2 ) + k * ′ 32 x 2 + l * 32 u 2 .
(31) In order to unequivocally determine u 1 , u 2 and u 3 , the following equation should have a unique solution

2 0 -l * 13 -l * 21 1 0 0 -l * 32 1 u 1 u 2 u 3 = 2k * ′ 1 x 1 + l * 1 r -(k * 1 -k * 13 ) ′ x 3 k * ′ 2 (x 2 -x 1 ) + k * ′ 21 x 1 k * ′ 3 (x 3 -x 2 ) + k * ′ 32 x 2 .
The determinant of the square matrix above is 2l * 13 l * 21 l * 32 = 1, so even in this case the ideal inputs u 1 , u 2 and u 3 are well defined. In addition, using similar calculations as in the previous case, it is possible to show that ė1 + ė13 = A m (e 1 + e 13 ) ė21 = A m e 21 , ė32 = A m e 32 .

Moving beyond the analysis with ideal input, we have that in the presence of parametric uncertainties, the following result, which extends the parametrization in (Baldi and Frasca (2018)) to general graphs, holds. Theorem 2. Under Assumptions 1 and 2, for any pair of agents (j, i), the dynamics of the synchronization error

N i=0 a ij ėji = A m N i=0 a ij e ji + b j • (33) N i=0 a ij k′ ji x i + k′ j N i=0 a ij e ji + N i=0 a ij lji u i
holds independently of the network connection.

Proof 2. The parametrization (33) turns out to be independent on the network connection thanks to the fact that all error dynamics can be homogenized to the reference model (A m , b m ) via appropriate control gains. Therefore, the dynamics can be summed (cf. ( 29) or (32)) even in the presence of cycles and undirected links.

Given the parametrization (33), one might be tempted to say that the algorithm in Theorem 1 can be used straightforwardly with the Lyapunov function ( 27). However, some attention must be paid when doing this: in fact, the actual input u j may be not well defined for all time instants on general graphs. To explain this point, let us collect all inputs in (25) on the left-hand side, leading to

U [u 1 • • • u N ] T = [β 1 • • • β N ]
T for an appropriate square matrix U depending on the estimates l ji . On directed acyclic graphs, U can be made upper triangular, with positive weights in its main diagonal, thus U is always invertible. On general graphs, the invertibility of U depends on l ji1 . Despite this difficult analytic aspect, the simulations in the next section show that the algorithm in Theorem 1 can handle networks beyond Assumption 3, and U turns out to be invertible at all time instants. Remark 3. The expression after (31) reveals that the agent are 'fictitiously' inverting the matrix U by only using neighbors' information. How such inversion is robust to delays is an open problem worth of future investigation.

ẋi = 0 1 a 1i a 2i x i + 0 b 1i u i (34)
with coefficients and initial conditions as in Table 1. The matrices are given in terms of ẋ0 = A m x 0 + b m r for the reference model and ẋi = A i x i + b i u i , i ∈ {1, . . . , N } for the other agents. The other design parameters are taken as:

γ k = 3, γ l = 0.3, Q = diag(1, 3
), and all estimated gains are initialized to 0. Also note that sgn(l * i ) = 1, ∀i. The simulations are carried out for a sinusoidal reference r of frequency 0.2 rad/s and amplitude 1. For the acyclic network, the resulting synchronization is shown in Figs. 4 and7. All states converge asymptotically to the state of the reference model.

Table 1. Coefficients and initial conditions of agents

a 1 a 2 b 1 x 0 agent #0 -0.5 -1 1 [1 -1] ′ agent #1 -1 2 1 [1 1] ′ agent #2 -0.75 2.5 0.5 [-1 -1] ′ agent #3 -1.25 2 1.25 [-1 0] ′ agent #4 -0.5 1 0.75 [0 1] ′ agent #5 -0.75 1 1.5 [1 0] ′ agent #6 -1.5 2.5 1 [-1 1] ′
In the directed cyclic graph, two cycles are present (2-3-4 and 4-5-6). Using the same parameters as in the previous simulations, the resulting synchronization is shown in Figs.

5 and 8. It can be seen that synchronization is slightly faster, at the price of a larger magnitude of the input. Finally, for the undirected graph the synchronization is shown in Figs. 6 and9). It is observed that having bidirectional connections does not necessarily help synchronization: synchronization is slower than in the previous cases. We studied synchronization of uncertain agents via a MRAC formulation with distributed input. We showed that the parametrization derived for the acyclic case (Baldi and Frasca (2018)) can be extended to more general graphs. Despite a suitable Lyapunov function exists and ideal inputs (with ideal gains) might be well defined, it is difficult to prove that the actual inputs (with estimated gains) are well defined for all time instant. Simulations suggest so. Future work will include considering unmatched uncertainties [START_REF] Lymperopoulos | Adaptive control of networked distributed systems with unknown interconnections[END_REF]) and switching topologies by using adaptive switching strategies [START_REF] Sang | Adaptive control of piecewise linear systems: the state tracking case[END_REF]; [START_REF] Yuan | Adaptive asymptotic tracking control of uncertain time-driven switched linear systems[END_REF]).

  is the input, and A 1 , A 2 , A 3 and b 1 , b 2 , b 3 are unknown matrices of appropriate dimensions, with possibly A 1 = A 2 = A 3 and b 1 = b 2 = b 3 . Time index will be omitted whenever obvious. Consider the reference model ẋm = A m x m + b m r (2) where x m ∈ R n is the state of the reference model, r ∈ R is its reference input, and A m and b m are known matrices of appropriate dimensions, with A m being Hurwitz so as to have bounded state trajectories x m .

  Fig. 3. Acyclic, cyclic and undirected networksThe simulations are carried out on a directed acyclic network, on a directed cyclic network, and on an undirected network, as shown in Fig.3. Node 1 acts as the leader node and the reference model is indicated as agent 0. The agents are second-order linear agents in the canonical form

Fig. 6 .Fig. 9 .

 69 Fig.6. Undirected network: state/input for reference model (dashed) and agents (solid)

  u 1 (12) where k 21 , k 2 , l 21 are the estimates of k * 21 , k * 2 , l * 21 respectively. The scalar l * 2 does not need to be estimated, only its sign is needed. The adaptation law in (12) is derived via the dynamics of the error e 21 = x 2x 1 ė21 = A m e 21 + b 2 ( k′ 21 x 1 + k′

						21 . By
	taking the Lyapunov function
	V 21 = e ′ 21 P e 21 +tr	k′ 21 k21 γ k |l * 2 |	+tr	k′ 2 k2 γ k |l * 2 |	+	l2 21 γ l |l * 2 |	(14)

2 e 21 + l21 u 1 ) (13) with k21 = k 21k * 21 , k2 = k 2k * 2 , l21 = l 21l *

  )γ k b ′ m P (e 31 + e 32 )(e 31 + e 32 ) ′ l31 = -sgn(l * 3 )γ l b ′ m P (e 31 + e 32 )u 1 l32 = -sgn(l * 3 )γ l b ′ m P (e 31 + e 32 )u 2 (18) where k 31 , k 32 , k 3 , l 31 , l 32 are the estimates of k * 31 , k * 32 , k * 3 , l * 31 , l * 32 respectively. We derive the adaptation law in (18) via the dynamics of the error e 321 = e 31 + e 32 and the Lyapunov function V 321 = e ′ 321 P e 321 + tr

	k′ 31 = -sgn(l * 3 )γ k b ′ m P (e 31 + e 32 )x ′ 1
	k′ 32 = -sgn(l * 3 )γ k b ′ m P (e 31 + e 32 )x ′ 2
	k′ 3 = -sgn(l * 3
				e 31 (t) + e 32 (t)
				2
	+l 31 (t)	u 1 (t) 2	+ l 32 (t)	u 2 (t) 2	(17)
	and the adaptive laws

A companion paper(Baldi et al. (2018)) shows that appropriate parameter projection can guarantee invertibility of U .