Physiologically structured cell population dynamic models with applications to combined drug delivery optimisation in oncology

.


Introduction
Optimising drug delivery in the general circulation targeted towards cancer cell populations, but inevitably reaching also proliferating healthy cell populations imposes to design optimised drug infusion algorithms in a dynamic way, i.e., controlling the growth of both populations simultaneously by the action of the drugs in use, wanted for cancer cells, and unwanted for toxic side effects on healthy cells.
Towards this goal, we design models and methods, with optional representation of circadian clock control on proliferation in both populations, according to three axes [17,18]: a) representing the oncologist's main weapons, drugs, and their fates in the organism by molecular-based pharmacokinetic-pharmacodynamic equations; b) representing the cell populations under attack by drugs, and their proliferation dynamics, including in the models molecular and functional targets for the drugs at stake, by physiologically structured equations; c) using numerical algorithms, optimising drug delivery under different constraints at the whole organism level, representing impacts of multiple drugs with different targets on cell populations.
In the present study, two molecular pharmacological ODE models, one for oxaliplatin, and one for 5-Fluorouracil, have been designed, using law of mass action and enzyme kinetics, to describe the fate of these two cytotoxic drugs in the organism.An age-structured PDE cell population model has been designed with drug control targets to represent the effects of oxaliplatin and 5-Fluorouracil on the cell cycle in proliferating cell populations.The models for proliferating cell population dynamics involve possible physiological fixed (i.e., out of reach of therapeutic influence, but of potential use by adapting drug delivery to it) circadian clock control, and varying drug control to be optimised, connected or not with pharmacological models.c EDP Sciences, 2016 Concentrations of drugs, represented by outputs of ODEs, are assumed to be homogeneous in the cell populations under attack by cytotoxic drugs.The possibility to describe the effects of other drugs, cytostatic (including in this category anti-angiogenic drugs, considered as acting on the G 1 phase, choking its entries and slowing it down), is also presented, but not put in pharmacokinetic equations and actual simulations in this study, that is focused on the combination of 5-FU and oxaliplatin, a classic therapeutic association in the treatment of colorectal cancer.
We then set conditions to numerically solve drug delivery optimisation problems (maximisation of cancer cell kill under the constraint of preserving healthy cells over a tolerability threshold) by considering a trade-off between therapeutic and toxic effects.The observed effects on proliferation are growth exponents, i.e., first eigenvalues of the linear PDE systems, in the two populations, healthy and cancer.The solutions to an optimisation problem taking into account circadian clock control are presented as best delivery time schedules for the two drugs used in combined treatments, to be implemented in programmable delivery pumps in the clinic.

Molecular pharmacokinetics-pharmacodynamics
We represent drug fate in the organism by physiologically based pharmacokinetic-pharmacodynamic (PK-PD) ODEs for drug concentrations.Drug effects, outputs of the PK-PD system, will later be considered as control inputs in a physiologically structured cell population dynamic model representing cancer or healthy tissue proliferation with prescribed targets (see Section 3).In this Section, we consider two cytotoxic drugs, oxaliplatin and 5-fluorouracil, that act on cells by killing them unless they undergo repair, and a cytostatic drug, assumed to exert a blocking action at phase transitions of the cell cycle without inducing cell death.
Drugs are described by concentrations, from their infusion in the circulation, which is represented by infusion flows assumed to be externally controlled by programmable pumps, until their action at the cell level.Whenever relevant, physiological control of biological mechanisms by circadian clocks is represented here by cosine-like curves -in the absence of better identified physiological gating functions -describing gate opening at transitions between cell cycle phases (0: gate closed; between 0 and 1: gate opening).

Oxaliplatin pharmacodynamics: impact on the cell cycle
Oxaliplatin is an anticancer drug that has been in use for more than 12 years now in the treatment of colorectal cancer.It exerts its action on cells by creating irreversible oxaliplatin-DNA adducts that subsequently yield doublestranded breaks in the DNA.Although it has been shown to induce apoptosis [31] in a manner that does not show major dependency on cell cycle phases (with the restriction that of course the DNA is always more exposed, and thus more sensitive to oxaliplatin insults, in S-phase), other studies have shown that its main effect on cancer cells was likely to be mostly active by arresting the cell division cycle at the G 2 /M, and to a lesser extent, at the G 1 /S transitions [81,82], effects on apoptosis occurring only after a prolonged exposure time.In the sequel, we shall assess these two main actions, both induced by DNA damage.The effect on death rates will be assumed to be directly proportional to the damage induced to the DNA, whereas the effect on phase transitions, physiologically due mainly to p53 (not represented as such in this study), will be represented by a variable delay due to cells under repair, that will be represented by a specific subpopulation appended to each one of the main phases G 1 and S G 2 .

Oxaliplatin PK-PD: from infusion until DNA damage
The PK-PD model for oxaliplatin runs as follows A simulation of the model behaviour under square wave drug infusion flow of oxaliplatin following a circadian rhythm is shown on Fig. 1.FIGURE 1. Behaviour of the two systems protecting the cells against oxaliplatin.Upper track: oscillations of plasma proteins in the presence of infused plasma oxaliplatin (4 consecutive days of chemotherapy followed by 10 days of recovery); lower track: oscillations of the concentrations of cellular reduced glutathione in presence of Pt 4+ ion, active intracellular metabolite of oxaliplatin.Decaying track: modelled apoptosis effect of oxaliplatin; in the absence of repair, the density of free DNA decreases in the mean, leading cells to death.The parameters used in these simulations (see Section 4.4 for details) have been chosen to produce expected effects, and have not been experimentally identified.
The features of the model are: irreversible binding of oxaliplatin to plasma proteins, synthesized in the liver; simultaneous attack in the intracellular medium of free (unbound) DNA, the drug target, and reduced glutathione, that thus acts as a competitive target of oxaliplatin, shielding the DNA from oxaliplatin; continuous synthesis of reduced glutathione; possible return to equilibrium value of DNA obtained by excision repair enzymes.Apart from blood pharmacokinetic constants that are easily accessible, other parameters have been evaluated so as to produce likely behaviour for the drug in tissues.
When circadian rhythms are considered relevant (i.e., among teams who practice anticancer circadian chronotherapy, cf.[33, 34, 47-50, 52, 53]), peak phases j PLP and j GSH for tuned values K 0 and G 0 are fixed according to laboratory or clinical published data [55,78] adapted to human organisms (when observations were performed on mice, it was assumed that 0 Hour After Light Onset (HALO) = 8 h, local time, and that Humans, contrary to mice, are diurnal animals).In the simulations presented, care was not taken to optimise drug delivery according to a clinically fixed circadian schedule, since optimisation will be performed freely later, using numerical algorithms and with no a priori knowledge on optimal delivery time schedules.Only the parameters of circadian physiological control functions for the cell division cycle and for the metabolism of the considered drugs will be fixed based on such laboratory or clinical data.Note that in this physiologically based model, circadian rhythms are not mandatory in the pharmacological model and can be neglected by setting d = e = 0 (and z = h = 0 in the FU model below) if they are considered as irrelevant.The idea to introduce them (optionally in our setting) in the pharmacological model comes from previous collaborative works with a team dealing with cancer chronotherapy, Francis Lévi's at Paul-Brousse Hospital in Villejuif, France.However, even though it is possible to neglect these chronopharmacological effects, it is essential to represent differences between healthy and tumour cell populations, and these are based in our model on a difference between physiologically fast and pathologically slow responses to circadian clock stimuli in the two populations, respectively healthy and cancer.This point will be discussed in the conclusion.

5-FU pharmacodynamics and the cell division cycle
The pharmacodynamics of 5-FU is complex and multiply-targeted, as shown, e.g., in [56], resulting in non viable modifications of the RNA and of the DNA, but most of all in thymidylate synthase inhibition.We will consider here only its effects on its main target, thymidylate synthase, an enzyme that is essential in DNA duplication and S-phase specific.In the same way as direct damage to the DNA was considered as proportional to added cell death, we will assume here that thymidylate synthase degradation is proportional to 5-FU-induced cell death.

5-FU PK-PD: from infusion to intracellular damage
As in the previous section (2.1) for oxaliplatin, the PK-PD of the cytotoxic drug 5-fluorouracil (5-FU) is represented by a system of ODEs, controlled by a continuous infusion flow j(t), with the help of a continuous flow k(t) of folinic acid, a.k.a.leucovorin, a natural compound of the folate family, that stabilises the complex formed by FdUMP with its target thymidylate synthase (TS) [56].Also for 5-FU have been experimentally evidenced circadian variations of the target TS and of the main enzyme responsible for drug catabolism, namely dihydropyrimydine dehydrogenase (DPD), such catabolism being processed mainly (80%) in the liver [56], which has been considered here as a filter on the general circulation.On entering into a cell, the active form of 5-FU, FdUMP, may be effluxed out of the cell by an ABC transporter [65,73], before reaching its target, a known mechanism of resistance to 5-FU [65].
, accounting for circadian tuning (if z 6 = 0, h 6 = 0) of equilibrium values for tissue thymidylate synthase (TS) activity S 0 and liver dihydropyrimidine dehydrogenase (DPD) maximal rate l DPD .The variables represented are concentrations, either in blood, or in tissues, for: P: plasma 5-FU, with source term a continuous infusion flow j(t), detoxicated at 80% in the liver by DPD U: intracellular FdUMP (main intracellular active metabolite of 5-FU) Q: plasma leucovorin, also known as folinic acid, with source term a continuous infusion flow k(t) L: intracellular methylene tetrahydrofolate, active metabolite of leucovorin N: nuclear activation factor of the ABC transporter MRP8 (ABCC11), assumed to be sensitive to damage to the FdUMP target thymidilate synthase (TS) A: ABC transporter activity (MRP8, or ABCC11, for 5-FU) S: free thymidylate synthase (TS), with equilibrium value S 0 (including circadian modulation if z 6 = 0) B: reversibly FdUMP-bound TS (binary complex U⇠T S) T : irreversibly FdUMP-bound TS (stable ternary complex U⇠L⇠T S) (B + T )/S: ratio of bound to unbound TS The main features of this molecular PK-PD model are: active degradation of plasma 5-FU by dihydropyrimydine dehydrogenase (DPD) in the liver, considered as a filter placed in the general circulation; active cellular intake in the target tissue, together with active efflux by ABC transporter, of the active metabolite FdUMP of 5-FU; stimulation of ABC transporter activity by nuclear factor, sensor of damage to drug target TS; reversible binding of intracellular metabolite FdUMP to its target, forming a binary complex; irreversible binding of the same by fixation of methylene tetrahydrofolate to the binary complex, yielding a stable ternary complex.A variant of this model, where the nuclear factor inducing cell efflux was assumed to be directly triggered by the drug, which we think physiologically less likely than by its intracellular effects on its target TS, has been published in [52].When circadian rhythms are considered relevant for PK-PD, peak phases j DPD and j T S for tuned values l DPD and S 0 are fixed according to laboratory or clinical published data [71,83], with adaptation from mice to Humans, as mentioned in the previous Section.
A simulation of the model behaviour under square wave drug infusion flows of 5-FU and folinic acid following a circadian rhythm is shown on Fig. 2.

Mutual exclusion from simultaneous infusions
Due to the chemical forms, always prepared for intravenous infusions (oxaliplatin is prepared as an acid solution whereas 5-FU is usually presented in an alkaline form, hence a risk of precipitate in the infusion line), in which the FIGURE 2. Behaviour of the variables (concentrations in molecules) under a joint infusion in the general circulation of 5-FU (P) and folinic acid (a.k. a. leucovorin, Q), mimicking a repeated sequence of 4 consecutive days of chemotherapy followed by 10 days of recovery.The respective intracellular metabolites, U for FdUMP and L for methylene tetrahydrofolate (MTHF), act synergistically, binding free thymidylate synthase, abbreviated in the sequel as TS (S), to form firstly a reversible binary complex FdUMP⇠T S (B) and secondly an irreversible ternary complex FdUMP⇠MT HF⇠T S (T ).The ABC transporter MRP8 activity (A) is stimulated by a nuclear factor (N), itself triggered by bound TS.The last track shows the ratio of bound-to-unbound TS, a measure of damages to TS.The parameters used in these simulations (see Section 4.4 for details) have been chosen to produce expected effects, and have not been experimentally identified.
two drugs are available in the clinic, it is not advisable to deliver them simultaneously.This has led us to add mutual exclusion terms for their infusion flows in the combined delivery to be optimised, as will be mentioned in Section 4. Rather than adding constraints on drug delivery schedules that are actually in use in the clic, we preferred, for computing convenience, to represent the risk of precipitate formation by an added artificial "precipitate" variable in the plasma, with supposed catastrophic consequences for toxicity to a healthy cell population: resulting in a term d prec .z(t) to be added to the irrecoverable death rates d i (t) and d kti (t) in the cell population dynamic model (3.3) in Section 3.2, see below.This precipitation risk is systematically taken into account and avoided at the infusion line level in chronotherapy delivery schedules [51,64] and in the FOLFOX protocol [25].

Cytostatic drugs: control on cell cycle phase transitions
As regards cytostatic drugs, that have been named so because they slow down tissue proliferation but, contrasting with cytotoxic drugs, do not -at least at low or medium doses -kill cells by damaging them and thus do not need trigger repair mechanisms, we can represent their action as firstly proposed in [11] by inhibiting cell cycle phase transition rates at G 1 /S and G 2 /M checkpoints.This representation is in particular adapted to cyclin-dependent kinase inhibitors (CDKIs) that control these checkpoints.It may also be used to represent the action of growth factor receptor antagonists, that may be thought of as decreasing the boundary terms at the beginning of each cell cycle phase, and also of anti-angiogenic drugs, that choke a solid tumour population, resulting in proliferating cell population models for the division cycle in decreased boundary terms at the beginning of phases G 1 and S.This point will be made more precise in Section 3. We have already mentioned that oxaliplatin, a certainly cytotoxic drug, can also be considered as having effects on phase transitions, as though it were additionally also a cytostatic drug.
A way to represent control on the cell cycle by cytostatic drugs, for instance tyrosine kinase inhibitors (TKIs), with introduction of a quiescent cell compartment, has been proposed in [37], using proliferative and quiescent cell compartments.The equations are written as where p(t, x) is the density of proliferating cells of age x at time t, and Q(t) is the density of quiescent (i.e., non-proliferating, out of the division cycle) cells at time t, and the drug target here is f , rate of escape at mitosis towards the siding phase G 0 (quiescent cell compartment), this f to be enhanced by a cytostatic drug.The model was identified on the human Non Small Cell Lung Cancer (NSCLC) cell line PC-9 submitted to the TKI drug erlotinib.
In fact, the equation for the quiescent phase is just a linear equation and, rewriting the system for the proliferating phase as we can see that it is nothing but the classical McKendrick transport equation ( [62], see also [11], recalled in Section (3.1)) for one proliferating phase with a balanced control between mitosis ((1 f )b ) and enhanced "cell disappearance" ( f b + µ), i.e., way out of the proliferating phase towards either quiescence or death.In other cell cycle models with phases G 1 , S, G 2 and M [11], this representation of a target for cytostatic drugs should be placed between phases G 1 and S, resulting in a modification of the McKendrick equation : with the adjunction of a quiescent phase G 0 represented by 8 < : which is thus fed only by cells escaping from G 1 instead of processing into S-phase.
In the sequel, we will focus only on the action of cytotoxic drugs oxaliplatin and 5-fluorouracil; nevertheless it is also obviously possible to use this local representation of the dynamics of cytostatic drugs to study and optimise combinations of both types of anticancer drugs, for instance oxaliplatin, 5-fluorouracil and cetuximab, as it is now frequently the case in the clinic of metastatic colorectal cancer [13,51].

Cell cycle model with drug damage and repair phases 3.1. Age-structured models of the cell cycle in cell populations
Various models of the cell division cycle at the single cell level exist, some of them with remarkable descriptions of their molecular mechanisms, allowing to study entrainment of the cell cycle by the circadian clock based on upto-date physiological knowledge [40][41][42]44], but they are not adapted to describe proliferation in cell populations, except by considering cellular automata, as in [1][2][3][4].To the aim of describing proliferation in tissues, we rather advocate age-structured partial differential equations (PDEs), that take into account the main source of variability in such proliferating cell populations, i.e., age in the division cycle.We recall here the McKendrick [62] (or Von Foerster-McKendrick) model for the cell division cycle in proliferating cell populations [11,[19][20][21][22][23][24]67], where n i (x,t) is the density of cells in cell cycle phase i of age x at time t: together with an initial condition (n i (t = 0,.))1iI .It may be shown [67] that for constant or time-periodic control on the d i and K i!i+1 , solutions to (3.1) satisfy n i (t, x) ⇠ C 0 N i (t, x)e lt , where C 0 is a real positive number and N i are eigenvectors of the eigenvalue problem The importance of the first eigenvalue l , which, provided that the system is growing, is positive and simple [67], and obviously governs (since the positive eigenvectors N i are bounded) the time-asymptotic behaviour of the cell populations n i (t, x) in each phase i, is thus fundamental, and it may be studied as a resulting observable, a function of the theoretical drug control conditions exerted onto the system.This has been done in [8][9][10][11].
In this system of transport equations, one for each phase i of the cell division cycle, the control on phase transition rates K i!i+1 (t, x)(i = 1, 2) may be decomposed in k i (x).y i (t), using an experimentally identified steplike function k i of age x [11] and a theoretical 24-hour periodic function y i of time t.The phase difference between y 1 and y 2 is set to a half-period (12 hours) because of the known phase opposition between the clock-controlled proteins Wee1 and p21 [43,61], which themselves control Cdk1 (that determines the G 2 /M transition) and Cdk2 (that determines the G 1 /S transition), respectively.Another (numerical) reason is that, as mentioned in [11], this phase difference has been found to maximise the growth coefficient, first eigenvalue of the periodically controlled McKendrick system; otherwise said, this value of a half-period for the phase difference between y 1 and y 2 yields the lowest population doubling time.Note that whereas in the PK-PD model circadian clocks are optional, in the cell cycle model they are fundamental to make the difference between healthy and cancer cell populations, as presented in Section 4.1.

A new cell cycle model with repair
In the sequel, we propose a more elaborate model, starting from the previous ones, for the cell population model under insult by drugs, taking into account various modes of action of cytotoxic drugs on the cell cycle, and knowing that the parameters of this model are many and will take many efforts to be experimentally identified.Nevertheless, we present it as a possible proof of concept for therapeutic optimisation procedures, based on representations that are partly molecular (for intracellular drug effects) and partly phenomenological (for the cell cycle), but always physiologically based, after our knowledge of the mechanisms that govern progression in the cell division cycle.
The proposed model, a variation of the previous McKendrick-like models for the division cycle in cell populations mentioned above, runs as follows and In this cell cycle model, 5 cell subpopulations are represented by their density variables: n i (i = 1, 2, 3) for agestructured phases G 1 , S G 2 and M, respectively, plus r i (i = 1, 2) for supplementary repair phases R 1 and R 2 (in which age is unchanged) corresponding to cells subject to DNA damage occurring in phases G 1 and S G 2 respectively.The durations of phases G 1 , S G 2 , R 1 and R 2 are not fixed (no upper bound in age, although in fact, they will never last infinitely long), whereas M-phase is assumed to have fixed duration (e.g., x M =1 hour) and to always (just by setting high the value of the constant M in the last equation, sending all cells to division at age x M in M-phase) lead in such fixed time each cell in mitosis into 2 intact daughter cells in G 1 .We assume that cells in M-phase are endowed with a supercoiled DNA, which make them unreachable by DNA-targeted drug insults, and we have not considered here the effects of M-phase-specific drugs (vinca alkaloids or taxanes); it is certainly possible to represent them by targets at the M/G 1 transition (i.e., by blocking mitosis), but thus far, the M-phase is free from pharmacodynamic effects in this model.
As regards the two cytotoxic drugs at stake in this study, apart from this last point about M-phase, oxaliplatin is a non phase-specific cytotoxic drug, whereas 5-FU is S-phase specific.Functions L i (t) quantify the effect of oxaliplatin (and of 5-FU for S G 2 ) on the cell subpopulation under attack in G 1 (i = 1) and S G 2 (i = 2) that can be coped with by repair, whereas d kti (t) functions stand for the cytotoxic effects of oxaliplatin that cannot be repaired (one can think of them as double-stranded breaks in the DNA [31]), be they occurring in G 1 , R 1 , S G 2 or R 2 .Note that we have not considered in our simulations (see Sections 3.2 and 4.4, where D 0 2 = 0) irrecoverable effects of 5-FU in our cell population model for the division cycle, assuming that its main target is apoptosis induced in S-phase, probably due to lack of thymine [36], but that these effects are amenable to repair.We do not consider in this model the intra-S checkpoint.Thus functions L i (t) and d kti (t) together represent the ultimate pharmacodynamics of cytotoxic drugs on their actual targets: the cell division cycle.
In each repair phase R i , fed from either G 1 or S G 2 by the corresponding "leak" function L i , damaged cells with density r i (t, x), keeping their age x unchanged, undergo a repair mechanism represented by a time constant e 1 i .Those cells that have been successfully repaired move back to the corresponding normal cycle phase, G 1 or S G 2 , and having become normal cells (n i ) again, are candidates to proceed to next phase at the corresponding checkpoint, which is controlled by functions k(x).yi (t) (multiplied or not by cytostatic drug pharmacodynamics (1 g ksi (t)), see below).At the same time, a more destructive cytotoxic death process represented by functions d kti (t) goes on, and is enhanced by more extended DNA damage, which is represented by a nonlinear effect as a function of DNA damage for functions d kti (t), that are tuned according to the cytotoxic drug at stake.We consider here that repair in phase R i consists only in a loss of time e 1 i , that does not depend on the drug dose ; however, the drug dose does influence the instantaneous rate L i (t) of cells sent to repair.A simplifying assumption made is that repair is always performed ad integrum, i.e., we have not considered here the possibility to represent errors in repair that may contribute, together with inefficient control of DNA at cell cycle checkpoints, to the known genomic diversity and instability, let alone aneuploidy, of cancer cell populations.Moreover, the checkpoints in this model represent cell cycle arrest by cumulative effects of oxaliplatin seen as sending cells to the repair phase, this drug having been reported to act more by arrest at transitions than by induction of apoptosis [81,82].
Normal transition rates (in the absence of damage to the cell), K 1 (t, x) at G 1 /S and K 2 (t, x) at G 2 /M are described in a multiplicative way by k i (x), the physiological transition rate in the absence of circadian control (determined by the distribution of cell cycle phase durations in the cell population, as shown in [11]), multiplied by y i (t), which represents circadian gating by p21 at G 1 /S and by Wee1 at G 2 /M [43,61].Thus physiologically, as already mentioned in Section (3.1),K i (t, x) = k i (x).y i (t) .
Additional cytostatic drugs, e.g.growth factor inhibitors, cyclin-dependent kinase inhibitors (CDKIs), or even antiangiogenic drugs, may be represented in this model for a cancer cell population by an inhibitory effect on G 1phase by choking boundary terms at cell cycle transitions M/G 1 or G 1 /S.Indeed, since these drugs are assumed to be non-cytotoxic, taking into account their effect should be done not on death terms, but only on the boundary terms, either n 1 (t, x = 0) or n 2 (t, x = 0), i.e., in the latter case on k 1 (x).y 1 (t) by a multiplication by (1 g ks1 (t)), where g ks1 (t) stands for the pharmacodynamics of such cytostatic drugs (of maximum value 1).This may be achieved by the adjunction of a quiescent phase, as in equation (2.3), that has been studied in detail with parameter identification on cell cultures in [37], using in our case the suggested modifications (2.5, 2.6) of the system.Note that another possible way of representing the action of cytostatics is by slowing down a non constant velocity factor dx dt of cell progression in G 1 introduced in the argument of the ∂ ∂ x operator, as proposed in [45].

Discretisation scheme for the cell population dynamics model
To perform numerical simulations, we discretised the PDE model (3.3).We adapted the discretisation scheme presented in [11], Section 3.3 by adding repair phases, transitions from the physiological phases towards these new phases and permanent transitions from these new phases towards the physiological phases.The number of cells in the discretised model is represented by n k, j i = n i (kDt, jD x) and the dynamics of this discretized variable is given for the time step k by n k+1 = M k n k , where n k is the concatenation according to ages j and phases i of the n k, j i s.The matrix M k is a block matrix, where each block represents transitions from one phase to the other.The size of each block is related to the duration of the corresponding phase.For instance, we fixed the duration of mitosis to one hour, so that we do not need to consider ages larger than one hour in M-phase.We give in Figure 4 the position of the nonzero elements for a given matrix M k .
The coefficients of the PDE model depend on the drug infusion schemes.The coefficients of the transition matrices M k will also depend on them accordingly.

Therapeutic optimisation 4.1. Optimal long term viable chemotherapy infusion schedules
Focusing on cancer chemotherapy, we propose to minimise the growth rate of the cancer cell population while maintaining the growth rate of the healthy cell population above a given toxicity threshold L (we chose in the ap- plication presented below L = 0.021h 1 ).Infusions here should be thought of as referring to the drugs Leucovorin, 5-FluoroUracil (5-FU) and oxaliplatin (lOHP) as explained in Section 2.
We consider two cell populations with their respective dynamics.We model both of them using the same agestructured cell population dynamics model (3.3) but with different parameters.In fact, we assume that the transition and death rates are the same when there is no circadian control, but that cancer cell populations show a looser response to physiological circadian control.More precisely, we model this phenomenon by different timedependencies for gating at transitions in (3.4) recalling from [11] that e is here a very small positive number (typically 10 10 ) put in the equations only to ensure y H i > 0 and y C i > 0, which may be shown as sufficient to imply irreducibility of the matrix M in numerical simulations (see below) and thus applicability of the Perron-Frobenius theory.
We assume that the drug has the same effect on both populations, i.e.L 1 (t), L 2 (t) and that death rates are the same.This modelling choice couples the behaviours of the populations through the drug infusions.
Then we represent the growth rates of both populations by the Floquet eigenvalue of the corresponding model (l C for cancer cells and l H for healthy cells).We shall search for an optimal drug infusion schedule among the set of measurable and bounded (classical assumptions in control theory [70] for instance) and T -periodic functions.
We obtain the following Floquet eigenvalue optimisation problem with constraints: g T -periodic, mesurable and bounded

24h-periodic drug infusions
The McKendrick population dynamics yields a cell population described by the number of cells n i (t, x), where i is the phase, t is the time and x is the age in the phase.As seen in Section 3.1, if the death and transition rates d i and K i!i+1 are T -periodic, then the solutions n i (t, x) to (3.1) are asymptotically -i.e., for large time t -equivalent (in a L 1 sense that can be made very precise [67]) to C 0 N i (t, x)e lt , where the eigenvector N i is T -periodic in the time variable.
Hence, searching for infusion strategies that optimise the Floquet eigenvalue of a biological system is a way to study its long term behaviour under drug infusion.Indeed, most of chemotherapy schedules are repeated until the patient responds to the treatment, and conversely stopped if the patient does not respond after a given number of chemotherapy courses.In this model, we take into account this periodic behaviour directly in the model and in our objective function.
Classically, chemotherapy schedules consist of several days of severe treatment followed by a longer duration of rest (2 days of treatment every other week for FOLFOX therapy [25], 4 days of treatment every other week or even 5 days of treatment every third week for chronotherapy [64]).We chose, instead of imposing a rest period, to add a toxicity constraint that imposes in a dynamic way that the growth rate of the healthy cell population always remains high enough for the patient to live with a reasonable quality of life, even though he may experience drug-induced diminished physiological functions.In order to be compatible with the periodicity of physiological circadian control, the periodicity of the treatment should be a multiple of 24h.Numerical tests that we have performed on our model show that a period of one day seems to be a good choice, because other locally optimal solutions that we found assuming longer periods finally proved to be all 24h-periodic.
But since we have coupled the population dynamics with the PK-PD model described in Section 2, we have to take care that periodic infusion of drugs in the same way results (asymptotically) in periodic transition and death rates in the population dynamics model.In the case of linear systems, the Floquet-Lyapunov theory (see [63] for instance) gives the following condition: if the eigenvalues of the so called monodromy matrix (also called Floquet multipliers) all have modulus smaller than 1, then the system is asymptotically stable, and a periodic input asymptotically results in a periodic output.
One can generalise this result to a T -periodic nonlinear system.For a given T -periodic control u let F be the map that to an initial state y 0 associates the state y T of the system at time T .If F is a contraction in the sense that there exists c < 1 such that for all initial states y 0 and z 0 , kF(y 0 ) F(z 0 )k = ky T z T k  cky 0 z 0 k, then a periodic input asymptotically results in a periodic output.We numerically checked this fact during computations.

Optimisation procedure
Our goal, formalised in (4.1), is to minimise the growth rate of the cancer cell population while maintaining the growth rate of the healthy cell population above a given toxicity threshold L .To obtain the numerical resolution of problem (4.1), we firstly discretise the PK-PD models (2.2) and (2.1) and the population dynamics model (3.3).This is the so-called direct method, which consists of a total discretisation of the control problem and then of solving the finite dimensional optimisation problem obtained.
In our setting, at each time, the control g l (t) is the flow of drug l (1  l  n drugs ) infused at time t.Here, n drugs = 3 for 5-FluoroUracil, leucovorin and oxaliplatin.The discretised control (g k l ) l,k will be the array of the infusion time step by time step and drug by drug.As we search for 24h-periodic controls with one drug, we only need to define g on one day, i.e. g 2 R n drugs ⇥N T , where N T = T /Dt is the number of time steps.
We need to bound the flow of each drug to avoid evident toxicity.We bounded the flow of 5-FU and leucovorin to ḡ1 = ḡ2 = 500 mg/m 2 /h and the flow of oxaliplatin to ḡ3 = 50 mg/m 2 /h.Hence, Given a discretised infusion strategy g, we build the matrices M C (g) and M H (g), that are the nonnegative matrices modelling the discretised dynamics of each cell population under drug infusion g when drug concentrations represented by the PK-PD models have reached their day-periodic asymptotic behaviour.
Thus, we study the optimisation of the growth rate in the discretised model, represented by the principal eigenvalue r of the nonnegative matrices M C for cancer cells and M H for healthy cells where M H = M N T H .. .M 1 H is given in Section 3.3 (see also [11], Section 3.3) and M C is defined accordingly.The discretised approximation of l C (g(•)) is then 1  T log r(M C (g)) , and 1 T log r(M H (g)) for l H (g(•)).We obtain the finite dimensional optimisation problem: We chose a discretisation step of 10 minutes because it may be considered as a lower limit to the half-life time of 5-FU and oxaliplatin in the plasma [12,26,68], which is most likely even lower than the half-time of its downstream molecular effects at the cell level, our concern here.The oldest ages represented in the discretisation scheme are 10 days for each phase except for mitosis (M-phase) where we choose 2h.
By the Perron-Frobenius theorem [6], we know that if M is nonnegative and irreducible, its principal eigenvalue r(M) is positive and is a simple eigenvalue.Moreover, the principal eigenvector is unique up to normalisation and can be chosen such that u(M) 0. In our setting (irreducibility of the matrix and the Perron-Frobenius theorem), one can naturally define a function r from the set of nonnegative and irreducible real matrices in R n into R + , that to a matrix associates its principal eigenvalue.General eigenvalue optimisation of non symmetric matrices is a difficult (non convex, non differentiable) problem: see [54] and [66] for two algorithms dealing with this problem.However, for nonnegative irreducible matrices, as the principal eigenvalue is simple, this implies that r is differentiable.Indeed, denoting by v and u the left and right eigenvectors of a matrix M associated with a simple eigenvalue r, the derivative of r at M can be written as [46]: Thus, as the objective function is differentiable, differentiable optimisation theory applies.We then get the complete gradient of the objective by the chain rule.We have a composite objective of the type where r is the principal eigenvalue, M H is the discretised population dynamics, L is the output of the discretised PK-PD and g is the discrete periodic control.To compute the gradient of the objective function, we thus compute the derivatives of all these functions, either in closed form (for log and M H ) or by an iterative algorithm (for r and L(g)) and we multiply them all up.Following [32], we set the stopping criterion of internal iterative algorithms according to the theory of consistent approximation [69].
To solve the non convex problem (4.2), we use the method of multipliers [7], which solves a sequence of non constrained optimisation problems whose solutions converge to a local optimum of the constrained problem (4.2).

Simulations
We coded the discretised versions of the PK-PD models (Section 2) and of the population dynamics model (Section 3.2), and the method of multipliers [7] in Scilab [76].
Locally optimal infusion strategy with a combination of leucovorin (dash-dotted line), 5-FU (dashed line) and oxaliplatin (solid line).These infusions are repeated every day in order to minimise the growth rate of the cancer cell population while maintaining the growth rate of the healthy cell population above the toxicity threshold of 0.021h 1 .
where j i is the probability density function of a Gamma law given by G is the gamma function, a 1 = 8.28, b 1 = 1.052h 1 , g 1 = 0h, a 2 = 3.42, b 2 = 1.47h 1 , g 2 = 6.75h (we substracted the fixed duration of 1 hour for the mitosis phase M from the duration of S G 2 M identified in [11], thus obtaining for j 2 an estimate of the probability density function of the duration of S G2, so that the index i = 2 actually stands for the G 2 /M transition).
At the end of the optimisation algorithm, we obtained a locally optimal solution for the discretised optimisation problem (4.2) and we reported it in Figure 6.
It consists of infusing the three drugs in the following order: leucovorin between 10 h 30 and 21 h 20 (stopping the infusion at 19 h in oder to avoid mixing medicines has a marginal impact on the performances), 5-FU between 19 and 20 h and finally oxaliplatin between 20 and 22 h at maximal flow.Then, these infusions are repeated every day until the situation justifies to stop the treatment.
Note that if we keep the circadian clock response differences between healthy and tumour populations, as illustrated on Figure 4.1, essential to define the behaviours of the two populations in our model, but forsake circadian control on the pharmacology of the drugs, we obtain a very close result, as illustrated on Figure 4.4.This point will be discussed in the conclusion.
Without drug infusion, the growth rate of cancer cells (0.0265h 1 ) is larger that the one of healthy cells (0.0234h 1 ), this difference being due to different responses to circadian clock control.This gives an evolution of the respective populations, cancer cells becoming more and more present: see Figure 8.By following the infusion strategy numerically determined by the optimisation algorithm, we obtained that the asymptotic growth rate in the healthy cell population was actually above the chosen toxicity threshold (l healthy L = 0.021h 1 ) and that the asymptotic growth rate of cancer cells was weakened to (0.0229h 1 ).Thus by this optimal infusion strategy, the asymptotic growth rate for cancer cells was reduced by 13.3 % while the growth rate for healthy cells was reduced by only 10.0 %, so that the toxicity constraint is satisfied.This gave us a description of the evolution of the respective populations, which is illustrated on Figure 9. FIGURE 7. The same optimal strategy as in Figure 6, here in the absence of circadian control on the molecular PK-PD model.Locally optimal infusion strategy with a combination of leucovorin (dash-dotted line), 5-FU (dashed line) and oxaliplatin (solid line).These infusions are repeated every day in order to minimise the growth rate of the cancer cell population while maintaining the growth rate of the healthy cell population above the toxicity threshold of 0.021h 1 .One sees that differences from Figure 6 are hardly perceptible, showing a limited impact of circadian chrono-PK-PD and conversely stressing the major importance of circadian physiological control on the cell cycle as regards differences between healthy and cancer cell populations in our model.
These infusions cause possibly recoverable damage to the cells (be they healthy or cancerous), sending them to repair, at a rate L 1 (t) in G 1 -phase and at a rate L 2 (t) in phase S G 2 (see Figure 10).
If we compare L 2 (t) (Figure 10) and the fraction of the cell populations in phase S G 2 (Figure 11), we can see that the best time to cause damage in phase S G 2 is when the proportion of cancer cells is significantly larger than the proportion of healthy cells in this phase, i.e., between 22 and 24 h.The strategy takes into account the delay between the infusion and the effective action of the drug at the cell level due to the pharmacokinetics and the pharmacodynamics (diffusion, metabolism and elimination) of the drugs.We may also remark that although the infusion strategy stops infusions for 13 hours per day, the drugs remain active all day long.
In our cell population dynamics model, the fraction of cancer cells in G 1 -phase is never significantly higher than the fraction of healthy cells in G 1 -phase, so that we cannot tune a significant difference between the two populations in G 1 , and L 1 (t) should remain small to limit the toxicity of the treatment.Hence, 5-FU, that causes little damage to cells in G 1 -phase, is the main drug in this infusion strategy.
Simulating the transition from the stationary state without drug to the stationary state with periodic drug infusion (Figure 12), we see that after a transition of around 10 days, the treatment performs as expected (l cancer = 0.0229h 1 and l healthy = 0.021h 1 ).As mentioned above, by this optimal infusion strategy, the asymptotic growth rate for cancer cells was reduced by 13.3 %, whereas the asymptotic growth rate for healthy cells was reduced by only 10.0 %, and the toxicity constraint l healthy L = 0.021h 1 was satisfied.Note also that the growth rate in the healthy cell population always remains positive, which means that the healthy cell population production never decreases.This does not mean that there is no negative drug effect on this healthy cell population: one should have in mind that we have represented here only the rate of production of healthy cells, e.g., in the gut or in the bone marrow, while for obvious homeostatic reasons, constant consumption (elimination in the intestinal lumen, or release in the general blood circulation) must be considered and should result here, together FIGURE 8. Evolution of the population of cancer (blue, above) and healthy (green, beneath) cells without drug infusion during 12 days.We can see that the populations have different exponential growth rates (l cancer = 0.0265h 1 and l healthy = 0.0234h 1 ).Cancer cells proliferate faster than healthy cells.with this decrease of production, in a decay in physiologically functional healthy cell population (in our example, mature villi enterocytes or mature blood cells).However, since the growth rate in the healthy cell population never becomes negative, the source of healthy cell production is preserved by the optimisation procedure.
The proposed infusion strategy combines the three drugs.This means that we obtain better results in terms of the trade-off between damage to cancer cells and toxicity to healthy cells by infusing smaller doses of each drug than a larger dose of one single dose.This confirms (once more) that the use of this combination is efficient for the fight against cancer.We shall remark however that the proposed strategy is different from the ones used today in the clinic.
Firstly, the optimisation algorithm yields a daily schedule consisting of infusions during 11 hours and rest periods of 13 hours whereas one usually infuses larger doses and then gives longer rest periods to the patient.As discussed in Subsection 4.2, the schedules that seem to take the best advantage from the influence of circadian clocks on the proliferation of cells have a period of one day.However in this work, we did not study the emergence of drug resistances among the cancer cell population, although resistance is one of the main causes of the failure of chemotherapies in the long term.One may wonder here whether resistances are more or less likely to appear when performing short, rather than long, rest periods.
Secondly, the locally optimal strategy returned by the optimisation algorithm suggests infusing oxaliplatin after 5-FU.In the FOLFOX schedule for instance [25], it is advised to infuse oxaliplatin before 5-FU.However, in [35], the authors show that infusing oxaliplatin before, during or after the exposure to 5-FU does not produce significant differences.The results of the optimisation procedure performed on our preliminary model (in which parameters have not been experimentally identified) suggest that infusing oxaliplatin after 5-FU may have some advantages.

Conclusion and future prospects
We have presented in this study what is to our best knowledge the first physiologically based model of action of a combination of cytotoxic drugs, commonly used in the clinic of colorectal cancer, on the division cycle -with consideration of the two main phase transition checkpoints G 1 /S and G 2 /M -in populations of cells, together FIGURE 9. Evolution of the population of cancer (blue, above) and healthy (green, beneath) cells with the drug infusion, starting at time 0, given by the algorithm.Healthy cells keep multiplying (l healthy = 0.021h 1 ) while the cancer cell population is weakened (l cancer = 0.0229h 1 ).FIGURE 12. Daily mean growth rates for cancer (solid line) and healthy cells (dashed line) when starting drug infusions at time 0. After a 10-day transitional phase, the biological system stabilises towards the expected asymptotic growth rate.with the use of numerical optimisation to meet the question of maximising tumour cell kill under the constraint of preserving healthy cell population from unwanted toxic side effects.
The molecular PK-PD models are based on physiological knowledge, can evolve with new findings in this area and their parameters -that are thus far only theoretical, tuned to produce observable effects on their molecular targets -should be identified, at least partly, by biological experiments led in collaboration with teams of pharmacologists.Obviously enough, these models are quite complex, and some of their features (e.g., ABC transporters for 5-FU) may be forsaken if they are not proven to be important, in order to simplify experimental identification of model parameters.So far, these PK-PD models should be only considered as providing biological bases for a proof of concept of therapeutic control optimisation on the cell population dynamic model for proliferation.
This cell population dynamic model is a new version of a system of McKendrick-like equations, the characteristics of which must be adapted to represent either cancer or healthy cell populations.In this setting, we have chosen, as we had done previously [9][10][11], to use the chronotherapeutic paradigm, in which healthy and cancer cell populations differ only by good or bad circadian clock control on cycle phase transitions.This should indicate new tracks to optimise cancer chronotherapeutic schedules, when the classic drugs oxaliplatin and 5-FU are combined with other drugs, in particular with irinotecan and cetuximab [51].
Also of note, since the model relies on physiological principles, it is possible, when chronotherapeutics has not proved useful (which may be due to persistence of good circadian clock control in tumours, or to bad and unrecoverable control by the central circadian pacemaker of the suprachiasmatic nuclei on the whole organism, impinging also healthy cell populations) to use other pathological differences in cancer cell populations, compared with healthy tissues.For instance, adding to this model a physiological representation of the control by the p53 protein of phase transitions, DNA repair or apoptosis launching, based on triggering of the ATM protein by DNA damage and subsequent launching of p53 oscillations (as presented in [27][28][29][30]60]) should prove useful.Indeed, it is known that p53 is mutated and inefficient in about 50% of solid tumours [77], and thus instead of circadian control, it is possible to use impaired -or not -p53 control to account for differences between the two cell populations.Another expected benefit of the adjunction of a DNA damage-ATM-p53 model should be to assess the respective roles in cell death and in cell cycle phase transition blockade of the drug oxaliplatin, assumed to act by enhancing cell death but also, and primarily, according to [81,82], by blocking G 1 /S and G 2 /M phase transitions.
Another valuable remark about the essential physiolological (or physiopathological) difference in responses to circadian clock stimuli between healthy an cancer cell populations is that it may be related to a poor quality of intercellular communications in the latter.Indeed, it has been proposed for a long time, and recently updated by the same author [79,80], that impaired gap junctions in cancer cells are a constantly observed phenomenon that may account for the specificity -in particular heterogeneity -of cancer cell populations.In our model, as illustrated on Figure 4.1, the extended support of the cancer gating function y C , with a long tail, compared to the much sharper healthy gating function y H , may perfectly be considered to account for such poor synchronisation between cells in the cancer cell population due to their impaired gap junctions, and possibly to other biophysical mechanisms such as impaired mitochondria, as mentioned in [14].In this respect, response to circadian clock stimuli at transitions between cell cycle phases is essential in our model, but circadian chronopharmacology of anticancer drugs is not.
Obviously also, to theoretically investigate optimisation of therapeutic procedures using new drugs that are known by their action on physiological targets to hit cancer cell, and not, or less, healthy cell populations, is another possible application of this way of representing pharmacodynamic effects of anticancer drugs.It requires a PK-PD model for each drug, but the target cell population dynamic model of proliferation can be based on the same principles, and the optimisation procedure can be strictly the same.Note also that we have chosen to optimise eigenvalues, since the population model is linear.If a linear model proves inadequate, it is also possible to dynamically optimise therapeutics using cell population numbers, as done in [5,16], using the same optimisation principles, but not Perron-Frobenius-Floquet theory.For theoretical methods in optimal control applied to cancer therapeutics in a more general setting, we refer to the book [75] and to its foregoer [74].
Finally, one must note that in the present representation, healthy and cancer cell populations are considered as evolving separately, without any population coupling effect between them, except by drugs.The question of evolution toward drug resistance, in particular, and optimisation of therapeutics to circumvent it (see a review in [14]), has led to rather different models of cell populations evolving according to a resistance phenotype, in the framework of cell Darwinism, advocated, e.g., in [38,39] and recently studied from a deterministic point of view by integro-differential equations in [15,[57][58][59]72].In these models, coupling between cell populations is present, representing competition between cells for space and nutrients by a nonlocal logistic term (in a Lotka-Volterra-like setting) in the global -cancer and healthy -cell population.These models are amenable to describe evolution towards drug resistance of a phenotype representing between-cell variability in the cancer cell population.
Whether it is possible and relevant address in the same model such phenotype evolution towards drug resistance and circadian modulation on proliferation -two phenomena that do not occur at the same time scale, and not necessarily independently, as the same disorganising mechanisms, present in cancer, that trigger heterogeneity and drug resistance in cancer cell populations (also reviewed in [14]) may also trigger the disruption of circadian clocks -or the cancer tissue response to circadian inputs -is still an open question.
, accounting for circadian tuning (if d 6 = 0, e 6 = 0) of equilibrium values for tissue reduced glutathione (GSH) G 0 and plasma protein (PLP) K 0 concentrations.The variables represented are concentrations, either in blood, or in tissues for: R : Plasma oxaliplatin, with source term a continuous infusion flow i(t) K : Plasma proteins, with equilibrium value K 0 (circadian-tuned if e 6 = 0) C : Active tissue oxaliplatin (Pt 4+ ion) F : Free DNA, with equilibrium value F 0 : the target of the drug, that establishes DNA adducts by disulphur bridges with it, later creating double strand breaks G : Reduced glutathione (GSH), tissue shield against DNA damage induced by oxaliplatin, with equilibrium value G 0 (circadian-tuned if d 6 = 0).

FIGURE 3 .
FIGURE 3. Sketched linear model of the cell division cycle in proliferating cell populations.Variables n i (t, x) are densities of cells at time t and age in phase x in subpopulations, i.e., cell cycle phases G 1 (i = 1), S G 2 (i = 2) and M (i = 3).Moreover, two additional subpopulations have been added, R 1 and R 2 , with age-independent variables r 1 (t) and r 2 (t), respectively, to describe the fate of those cells that have been hit by drug-induced DNA damage and are waiting to be repaired -or sent to apoptosis -before being able to join next phase at G 1 /S and G 2 /M checkpoints, respectively.

FIGURE 4 .
FIGURE 4. Pattern of the nonzero elements of the transition matrix of the discretised dynamics.Each block represents the transitions from one phase to the other.The block G 1 /G 1 is a subdiagonal matrix representing the fact that at each time step, the cells that have stayed in G 1 get older.The block R 1 /R 1 is a diagonal matrix that represents cell proliferation arrest until cells go back to G 1 -phase .

FIGURE 5 .
FIGURE 5. Gating at cell cycle phase transitions due to circadian clock control.Functions y H i (left, healthy) and y C i (right, cancerous) define the hours at which the cells of each population (healthy and cancerous) can change phase: when y(t) = 0, transition gates are closed.The dashdotted lines correspond to the transition from G 1 to S G 2 and the solid lines correspond to the transition from S G 2 to M. In the model presented here, healthy and cancerous cells differ only by their responses to circadian clock control, represented by functions y H i and y C i .

FIGURE 10 .
FIGURE 10.Transition rate to repair phases R 1 and R 2 under the locally optimal infusion strategy presented in Figure6.The L 1 (t) and L 2 (t) entries in the population dynamics model (3.3) represent the (recoverable) damage on cells induced by the drug infusions.

FIGURE 11 .
FIGURE 11.Fraction of the cancer cell population in phase S G 2 (solid line) among all cancer cells and fraction of healthy cell population in phase S G 2 (dotted line) among all healthy cells.
. For healthy cells, we set K H i (t, x) = k i (x).yH i (t) and for cancer cells, we setK C i (t, x) = k i (x).yC i (t)where y H i (t) and y C i (t) are given, similarly to the ones in [11], by y H 1 (t) = cos 2 (2p(t 10)/12)1l [7;13] (t) + e,