Elisabetta De Maria

François Fages

Aurélien Rizk

Sylvain Soliman

Design, Optimization and Predictions of a Coupled Model of the Cell Cycle, Circadian Clock, DNA Repair System, Irinotecan Metabolism and Exposure Control under Temporal Logic Constraints 1

Keywords: model coupling, temporal logic, model checking, constraint solving, parameter learning, cell cycle, DNA damage, irinotecan

In systems biology, the number of available models of cellular processes increases rapidly, but re-using models in different contexts or for different questions remains a challenging issue. In this paper, we study the coupling of different models playing a role in the mammalian cell cycle and in cancer therapies. We show how the formalization of experimental observations in temporal logic with numerical constraints can be used to compute the unknown coupling kinetics parameter values agreeing with experimental data. This constraint-based approach to computing with partial information is illustrated through the design of a complex model of the mammalian cell cycle, the circadian clock, the p53/Mdm2 DNA-damage repair system, the metabolism of irinotecan and the control of cell exposure to it. We discuss the use of this model for cancer chronotherapies and evaluate its predictive power with respect to circadian core gene knock-outs.

Introduction

In systems biology, the number of available models of cellular processes increases rapidly. To date, most of the effort has been devoted to building models and making them freely available, through the design of standard exchange formats, such as for instance the Systems Markup Language SBML [START_REF] Hucka | The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models[END_REF], the making of model repositories, such as for instance Biomodels 2 , the making of biological ontologies to establish the links between molecular synonyms, species, units, etc., and the development of modeling tools, such as Cell Designer, Biocham [START_REF] Calzone | BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge[END_REF], BioNetGen [START_REF] Michael | BioNetGen: software for rule-based modeling of signal transduction based on the interactions of molecular domains[END_REF], Pathway Logic [START_REF] Eker | Pathway logic: Symbolic analysis of biological signaling[END_REF], Bio-ambients [START_REF] Regev | Representation and simulation of biochemical processes using the pi-calculus process algebra[END_REF], etc. Despite these efforts however, re-using models in different contexts or for different questions remains a challenging issue. In practice, most of the models are developed, refined, simplified or coupled with respect to other models by hand with no direct support from the tools to re-use models in a systematic way using a specification of the global behavior of the system.

Coupling biological models is necessary to study how the building blocks interact together and make predictions on the global system's behavior. Model coupling is also a method to better understand and improve the composite models. The knowledge acquired from the global view provided by a coupled model can indeed lead to modify the single model components in order to satisfy some observed property of the global system. In particular, coupling models can help identifying lacks in the model components, like a missing node in a pathway for instance.

In this paper, we show how the formalization of experimental observations in temporal logic with numerical constraints can be used to automatically find parameter values for the coupling kinetics agreeing with experimental data. We illustrate this constraint-based approach to computing with partial information, through the coupling of existing biochemical models of the mammalian cell cycle, the circadian clock, the p53/Mdm2 DNA-damage repair system, and irinotecan metabolism. Finally, we discuss the predictive power of the obtained coupled model with respect to circadian core gene knock-outs.

Mammalian cell cycle

Irinotecan is an anti-carcinogenic inhibitor of topoisomerase-1 which started to be used in clinical treatments approximately twenty years ago [START_REF] Nitiss | Dna topoisomerase-targeting antitumor drugs can be studied in yeast[END_REF]. It shows significant efficacy against a variety of solid tumors, including lung, colorectal, and cervical cancers. Scientists are currently trying to optimize the irinotecan therapy in order to understand how to limit its toxicity on healthy cells and to increase its efficacy [START_REF] Altinok | A cell cycle automaton model for probing circadian patterns of anticancer drug delivery[END_REF]. In this context, it is crucial to comprehend how the administration of this medicament influences cellular proliferation. For this purpose, the observed effects of the circadian rhythm on the toxicity and efficacy of anti-tumor drugs should be taken into account. In fact, the effectiveness of anti-cancer drugs on a healthy as well as tumorous cells is dependent on the phase of the cell cycle in which those cell lie [START_REF] Altinok | A cell cycle automaton model for probing circadian patterns of anticancer drug delivery[END_REF]. Under the hypothesis that the cell cycle in healthy tissues is mainly entrained by the circadian clock, it is possible to reduce the toxicity on healthy cells by injecting antitumor drugs in precise periods of the circadian clock.

On the other hand, tumorous cells are either phase-shifted (slow-growing tumors) or not entrained any more (rapidly growing or advanced stage tumors). A rhythmic drug exposure can thus limit toxicity on healthy cells while maintaining efficacy on tumour cells.

In this paper, we develop a complex model of the mammalian cell cycle, circadian clock, p53/Mdm2 DNA repair system and ironotecan metabolism to investigate the influences of irinotecan on cell proliferation. There are in the literature many models of the mammalian cell cycle [START_REF] Novák | A model for restriction point control of the mammalian cell cycle[END_REF][START_REF] Gérard | Temporal self-organization of the cyclin/cdk network driving the mammalian cell cycle[END_REF] and of the circadian biochemical clock [START_REF] Leloup | Toward a detailed computational model for the mammalian circadian clock[END_REF][START_REF] Daniel | A detailed predictive model of the mammalian circadian clock[END_REF], a few ones of the cell's DNA-damage repair network [START_REF] Ciliberto | Steady states and oscillations in the p53/mdm2 network[END_REF][START_REF] Chickarmane | A model for p53 dynamics triggered by dna damage[END_REF], and recently some preliminary models of irinotecan intracellular pharmacodynamics [START_REF] Dimitrio | Irinotecan: Modelling intracellular pharmacokinetics and pharmacodynamics. m2 master thesis[END_REF][START_REF] Ballesta | A combined biological and mathematical study of the anticancer drug irinotecan molecular pharmacokineticspharmacodynamics and their control by the circadian clock[END_REF]. However these modules need to be composed in a coherent way to make meaningful predictions.

Modeling under temporal logic constraints

Our approach to modeling in systems biology consists in formalizing the relevant properties of the behavior of the global system in temporal logic, and in using model-checking, constraint solving and continuous optimization algorithms to compute unknown parameters and validate the model with respect to its temporal specification. This temporal logic based approach is at the heart of our modeling platform, the Biochemical Abstract Machine Biocham [START_REF] Calzone | BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge[END_REF][START_REF] Fages | BIOCHAM v2.8 user's manual[END_REF].

Model-checking is the process of algorithmically verifying whether a given state transition structure is a model for a given temporal logic formula [START_REF] Clarke | Model Checking[END_REF]. In the literature, there are now various applications of model-checking techniques to biology. In [START_REF] Chabrier | Symbolic model checking of biochemical networks[END_REF][START_REF] Eker | Pathway logic: Symbolic analysis of biological signaling[END_REF], temporal logic was first introduced as a query language for biochemical networks and for validating boolean models of biological processes. Some experimental results were obtained on a large scale with Kohn's map [START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF] of the mammalian cell cycle control [START_REF] Chabrier-Rivier | Modeling and querying biochemical interaction networks[END_REF] (800 reaction rules, 500 variables) using the symbolic model-checker NuSMV, and on a small ordinary differential equation (ODE) model using the constraintbased model checker DMC. This approach to verifying biological processes has pushed the development of model-checking techniques for quantitative properties, and continuous, stochastic or hybrid models.

For (non-linear) ODE models, numerical integration techniques provide numerical traces on which formulae of Linear Time Logic with numerical constraints over R, named LTL(R) , can also be evaluated by model-checking [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF]. Simpathica [START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF] and Biocham are two computational tools integrating such model-checkers for quantitative models. This approach has been further developed in Biocham by generalizing model-checking to a temporal logic constraint solving algorithm [START_REF] Fages | From model-checking to temporal logic constraint solving[END_REF], allowing for efficient kinetic parameter optimization [START_REF] Rizk | On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology[END_REF] and robustness analysis [START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF] w.r.t. quantitative temporal properties formalized in LTL(R) [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF].

Related work concerns stochastic models and parameter uncertainty studies. In [START_REF] Heath | Probabilistic model checking of complex biological pathways[END_REF], Heath et al. apply the probabilistic model-checker PRISM to the study of a complex biological system, namely, the Fibroblast Growth Factor (FGF) signalling pathway. In [START_REF] Clarke | Statistical model checking in biolab: Applications to the automated analysis of t-cell receptor signaling pathway[END_REF], Clarke et al. apply statistical model-checking on a stochastic model of a T-cell receptor. In [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] Batt et al. develop a modeling framework based on differential equations to analyze genetic regulatory networks with parameter uncertainty. The values of uncertain parameters are given in terms of intervals and dynamical properties of the networks are expressed in temporal logic. Model-checking techniques are then exploited to prove that, for every possible parameter value, the modeled systems satisfy the expected properties and to find valid subsets of a given set of parameter values (such an approach is exploited in RoVerGeNe, a tool for robust verification of gene networks). In [START_REF] Piazza | Algorithmic algebraic model checking i: Challenges from systems biology[END_REF], Piazza et al. propose semialgebraic hybrid systems as a natural framework for modeling biochemical networks, taking advantage of the decidability of the model-checking problem for Timed Computation Tree Logic.

In this paper, we focus on the use of LTL(R) temporal constraints for integrating biochemical models. In order to compose the different modules, we assume a finite set of hypotheses concerning the structure of the links. The unknown kinetic parameter values are then computed by solving the temporal logic constraints using an evolutionary continuous optimization algorithm [START_REF] Rizk | On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology[END_REF] in order to make the model components interact in a proper way. For this, the biological properties of the global system are formalized as LTL(R) constraints, and solved so that the expected properties are automatically satisfied by the coupled model.

Organization of the paper

The paper is organized as follows. In Section 2 we introduce the temporal logic with numerical constraints LTL(R) used to specify relevant properties of both the composite models and the coupled model. Section 3 describes the elementary cell processes considered, their models taken from the litterature and their specification in LTL(R) . Section 4 presents the coupling of these elementary models and the specification of the global properties of the system in LTL(R) . Section 5 gives some performance figures for the evaluation of the parameter optimization method and the inference of parameter values . Then Section 6 shows how the coupled model and the parameter search method can be used to derive an optimal control model for maximizing the volume of irinotecan under non-toxicity constraint in synchronized (healthy) cells. Finally, in Section 7 we illustrate the predictive power of our coupled model by investigating the effects of clock genes knock outs on the cell cycle in silico and comparing the results with the literature. All the models and the temporal logic formulae used are available in (SBML compatible) Biocham format at3 .

Preliminaries on rule-based modeling and LTL(R) temporal logic specifications

The Systems Biology Markup Language (SBML) is a widely used rulebased formalism to describe systems of biochemical reactions. SBML is a useful format for exchanging models between modelers, and has been adopted for large repositories of models, such as for instance biomodels.net.

The rule-based language of Biocham for describing reaction models is compatible with SBML. Biocham adds a specification language based on temporal logic for formalizing the global properties of the system observed in biological experiments, under various conditions or gene mutations. Having formal languages not only for describing biochemical reaction models, but also for specifying their behavior, opens a whole avenue of research for designing automated reasoning tools to help the modeler [START_REF] Fages | Temporal logic constraints in the biochemical abstract machine BIOCHAM (invited talk)[END_REF].

The biological properties of quantitative models can be formalized in Biocham by formulae of the Linear Time Logic with numerical constraints over the reals LTL(R) [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF][START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF][START_REF] Rizk | On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology[END_REF]. LTL(R) formulae are formed over firstorder atomic formulae with equality, inequality and arithmetic operators ranging over real values of concentrations and of their derivatives, using the logical connectives and the usual temporal operators of LTL(R) : in particular operator G for "always in the future", F for "sometimes in the future", the next time operator X, and the binary operator until U.

For instance, F([A] > 10) expresses that the concentration of A eventually gets above the threshold value 10 and G([A] + [B] < [C]) states that the concentration of C is always greater than the sum of the concentrations of A and B. Oscillation properties, abbreviated as oscil(M, K), are defined as

The abbreviated formula oscil(M, K, V ) adds the constraint that the maximum concentration of M must be above the threshold V in at least K oscillations while period(M, P ) states that M oscillates at least 3 times and has a period P for the last three oscillations. It is worth noting that this expression of oscillations in temporal logic does not impose us to fix the phase and period of oscillations as in curve fitting.

LTL(R) formulae are interpreted in linear state transition structures which represent either an experimental data time series or a simulation trace, both completed with loops on terminal states. Given the ODE corresponding to a reaction model, under the hypothesis that the initial state is completely defined, a discrete simulation trace can be obtained by means of a numerical integration method (namely Rosenbrock method for stiff systems). Since constraints refer not only to concentrations, but also to their derivatives, traces of the form (< t 0 , x 0 , dx 0 /dt >, < t 1 , x 1 , dx 1 /dt >, . . . ) are considered, where at each time point t i , the trace associates the concentration values x i to the variables, and the values of their first derivatives dx i /dt. It is worth noting that in adaptive step size integration methods of ODE systems, the step size t i+1 -t i is not constant and is determined through an estimation of the error made by the discretization. The notion of next state refers to the state of the following time point in a discretized trace, and thus does not necessarily imply a real time neighborhood. The rationale is that the numerical trace contains enough relevant points, and in particular those where the derivatives change abruptly, to correctly evaluate temporal logic formulae.

Beyond verifying whether an LTL(R) formula is satisfied in a numerical trace (model-checking), an original algorithm for solving LTL(R) constraints [START_REF] Fages | From model-checking to temporal logic constraint solving[END_REF][START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF] has been introduced to compute a continuous satisfaction degree in [0, 1] for LTL(R) formulae [START_REF] Rizk | On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology[END_REF], opening up the field of model-checking to optimization. This is implemented in Biocham using an evolutionnary continuous optimization algorithm [] for optimizing parameter values with respect to LTL(R) properties.

Elementary Cell Process Models and Temporal Specifications

In this section we introduce the biological processes we deal with, giving temporal logic formulae to specify the behaviour of each of them. Each property is expressed first in natural language, then formalized in LTL(R) .

Mammalian Cell Cycle Control

Cells reproduce by duplicating their contents and then dividing in two. To produce a pair of genetically identical daughter cells, the DNA has to be faithfully replicated, and the replicated chromosomes have to be segregated into two separate cells. The duration of the cell cycle varies greatly from one cell type to another; in many mammalian cells it lasts about 24 hours. The cycle is traditionally divided into the following four distinct phases [START_REF] Alberts | Molecular Biology of the Cell[END_REF]: the G1-phase, that is the temporal gap between the completion of mitosis and the beginning of DNA synthesis, the S-phase (synthesis), that is the period of DNA replication, the G2-phase, that is the temporal gap between the end of DNA synthesis and the beginning of mitosis, and the M-phase (mitosis), when replicated DNA molecules are finally separated in two daughter cells.

The cell cycle is regulated by different checkpoints, that are moments when the cell progression is stopped to verify the state of the cell and, if needed, to repair it before damaged DNA is transmitted to progeny cells. DNA damaging agents trigger checkpoints that produce arrest in G1 and G2 stages of the cell cycle. Cells can also arrest in S, which amounts to a prolonged S phase with slowed DNA synthesis. Arrest in G1 allows repair before DNA replication, whereas arrest in G2 allows repair before chromosome separation in mitosis.

The proper alternation between synthesis and mitosis is coordinated by a complicated network that regulates the activity of a family of key proteins. These proteins are composed of two subunits: a regulatory subunit, a cyclin, and a catalytic subunit, the cyclin-dependent kinase, cdk for short. A cdk has to associate with a cyclin partner to form a dimer and has to be appropriately phosphorylated in order to be active. The progression through cell cycle is orchestrated by the rise and fall of the Cdk/cyclin dimers which are characteristic of each phase.

In this work we refer to the model of mammalian cell division proposed by Novák and Tyson in [START_REF] Novák | A model for restriction point control of the mammalian cell cycle[END_REF] and extended by Zámborszky et al. in [START_REF] Zámborszky | Computational analysis of mammalian cell division gated by a circadian clock: Quantized cell cycles and cell size[END_REF] to include the regulatory activity of Wee1, a kinase that delays or prevents mitosis by phosphorylation of the Cdk1/CyclinB complex. The extended model comprises 22 differential equations and 4 steady-state relations. This is the specification for a 100-hours simulation of the model: F cell : CycA is greater than 2 in at least 4 oscillations and CycB is greater than 3.5 in at least 4 oscillations and CycD is greater than 0.4 in at least 4 oscillations and CycE is greater than 1 in at least 4 oscillations.

LTL(R) : oscil([CycA], 4, 2) ∧ oscil([CycB], 4, 3.5) ∧ oscil([CycD], 4, 0.4) ∧ oscil([CycE], 4, 1).

Mammalian Circadian Clock

In many living organisms, the activity of some genes and proteins spontaneously display sustained oscillations with a period close to 24 hours. A biochemical clock present in each cell is responsible for maintaining these oscillations at this period. In mammalian cells, two major proteins are transcribed by clock genes in a circadian manner, CLOCK and BMAL1, which in turn bind to form a heterodimer responsible for the transcription of PER (Period) and CRY (Cryptochrome). The two newly-formed proteins then bind as soon as the activity of the complex reaches a threshold. PER/CRY associates with the complex CLOCK/BMAL1 to inhibit its activity and therefore the transcription of the two proteins PER and CRY. This negative feedback loop gives rise to sustained oscillations.

The adaptation of biological organisms to their periodically varying environment is mediated through the entrainment of circadian rhythms by light-dark (LD) cycles. Light can entrain circadian rhythms by inducing the expression of the PER gene.

The model of the circadian clock considered in this work is the one proposed by Leloup and Goldbeter in [START_REF] Leloup | Toward a detailed computational model for the mammalian circadian clock[END_REF], that consists of 19 differential equations incorporating the regulatory effects exerted on gene expression by the PER, CRY, BMAL1, CLOCK, and REV-ERBα proteins, as well as posttranslational regulation on these proteins by reversible phosphorylation, and light-induced PER expression.

The cyclic behaviour of the main compounds of the system is specified by the following formula: F clock : mPER, mCRY, mBmal1, and mREVERB oscillate with a period equal to 24 (in the last three oscillations). LTL(R) : period(mP ER, 24) ∧ period(mCRY, 24) ∧ period(mBmal1, 24) ∧

period(mREV ERB, 24).

Some recent researches showed the existence of biochemical links between the circadian and the cell cycle. In particular, Matsuo et al. [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF] proved that a cell cycle regulator, Wee1, is directly regulated by clock components.

P53/Mdm2 DNA-damage Repair System

The third model is devoted to the description of protein p53, a tumor suppressor protein which is activated in reply to DNA damage. P53 has the capability to arrest the cell cycle in the different phases and to lead to apoptosis, i.e. cell death. P53 can be activated in many ways, in particular in response to DNA damage.

In normal conditions, the concentration of p53 in the nucleus of a cell is feeble: its level is controlled by another protein, Mdm2. These two proteins present a loop of negative regulation. In fact, p53 activates the transcription of Mdm2 while the latter accelerates the degradation of the former. DNA damage increases the degradation rate of Mdm2 so that the control of this protein on p53 becomes weaker and p53 can exercise its functions. This protein is responsible for the activation of many mechanisms: in an indirect way, it stops the DNA synthesis process, it activates the production of proteins charged with DNA reparation, and can lead to apoptosis.

When DNA is damaged, Mdm2 looses its influence on p53 and one can observe oscillations of p53 and Mdm2 concentrations. The response to a stronger damage is a higher number of oscillations. Oscillations have a very regular period. In literature, several models have been proposed to model the oscillatory behaviour of proteins p53 and Mdm2, most notably the ones proposed by Chickermane et al. [START_REF] Chickarmane | A model for p53 dynamics triggered by dna damage[END_REF], by Ciliberto et al. [START_REF] Ciliberto | Steady states and oscillations in the p53/mdm2 network[END_REF], and by Geva-Zatorsky et al. [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF]. In this work we build upon the one described in [START_REF] Ciliberto | Steady states and oscillations in the p53/mdm2 network[END_REF], that consists of 6 differential equations.

The following three properties concern the behaviour of proteins p53 and Mdm2. F 1 p53 : In case there is no DNA damage, p53 and Mdm2 are constant functions.

LTL(R) :G([DN Adam] = 0) → G(d([p53])/dt = 0 ∧ d([M dm2 :: n])/dt = 0).

Irinotecan Metabolism

Camptothecins are substances that can be extracted from the Chinese tree "Camptotheca acuminata Decne" and are mainly used for the treatment of digestive cancers. Their anticancerogenic properties have been discovered at the end of the Fifties but the first clinical tests have been interrupted owing to heavy effects due to the toxicity of the substances. In the Eighties researchers discovered that camptothecins are inhibitors of topoisomerase-1 (Top1 for short), an essential enzyme for DNA synthesis. Afterwards, they started to focus on some semi-synthetic derivative of water-soluble camptothecins, such as irinotecan and topotecan. Irinotecan is pro-medicine and must be transformed in its active metabolite, SN38, to be effectively cytotoxic. In fact the anticancerogenic activity of irinotecan (CPT11) is approximately 100 times less effective than the one of SN38. The activation is due to carboxylesterase, an enzyme mainly located in the liver, in the intestine, and in the tumoral tissues. SN38 is then detoxified through glucorono-conjugation: this realizes uridine diphosphate glucoronosyl transferase 1A1. Mechanisms through which irinotecan damages the cell are very complex and have not been completely explained yet. It is sure that DNA lesions appear after the inhibition of Top1 by SN38. Top1 is a protein which is present in all living organisms and which checks DNA replication and transcription. It intervenes to modify the DNA winding degree, acting on one strand. More precisely, Top1 links itself to the extremity 3' of DNA forming a transitory cleavage complex and cuts a DNA strand, that in such a way is able to unroll. Then such a complex dissociates and a new ligature comes up. In normal conditions, the connection process is favored with respect to the cleavage one. The target of irinotecan, and above all of its active metabolite SN38, is the complex Top1-DNA. SN38 links to the complex through a covalent bond, preventing in such a way from the ligature of the DNA strand. As clearly written in the title of [START_REF] Pommier | Camptothecins and topoisomerase i: A foot in the door. targeting the genome beyond topoisomerase i with camptothecins and novel anticancer drugs: Importance of dna replication, repair and cell cycle checkpoints[END_REF], SN38 acts like a "foot in the door": it keeps opened the DNA strand to which Top1 is linked as to prevent a door from closing. These complexes are still reversible and do not cause DNA lesions. However, they favor them: some lesions can rise as a consequence of the possible collisions with the transcription complexes or with the replication fork. This induces the arrest of the cell cycle. In this case we speak of irreversible complexes. Lesions due to the inhibition of Top1 are therefore consecutive to the stages of the cell metabolism. It means that irinotecan injections must be repeated and abundant in order to be effective. Besides irinotecan is more effective during the DNA replication phase [START_REF] Ohdo | Cell cycledependent chronotoxicity of irinotecan hydrochloride in mice[END_REF][START_REF] Zhou | Transcriptional regulation of mitotic genes by camptothecin-induced dna damage : Microarray analysis of doseand time-dependent effects[END_REF]. Fur-thermore, the inhibition of the DNA synthesis takes rapidly place (in a few minutes) and lasts several hours.

Defence answers of cells subjected to irinotecan injections are multiple and vary according to the drug dose. The administration of a very light dose suffices to slow down the S phase of the cell cycle and to delay the G2-M transition. If the dose is more substantial, the lag time in the S phase is much more significant and the cell cycle arrest in the G2-M transition can last more than sixty hours or even be permanent. In this latter case, some genes responsible for the cell cycle arrest (as an example, p21) and involved in the aptototic pathway are over-expressed. These genes are activated by p53, and this suggests the intervention of the protein in reply to a DNA damage due to the dissociation of Top1 from DNA [START_REF] Zhou | Transcriptional regulation of mitotic genes by camptothecin-induced dna damage : Microarray analysis of doseand time-dependent effects[END_REF].

In this work we refer to a pharmacokinetics/pharmacodinamics (PK/PD) model of irinotecan developed by Dimitrio [START_REF] Dimitrio | Irinotecan: Modelling intracellular pharmacokinetics and pharmacodynamics. m2 master thesis[END_REF] and currently further elaborated by Ballesta [START_REF] Ballesta | A combined biological and mathematical study of the anticancer drug irinotecan molecular pharmacokineticspharmacodynamics and their control by the circadian clock[END_REF], that takes aim at representing the action of the drug on the body (pharmacodinamic) and the action of the body on the drug (pharmacokinetic), and thus the drug metabolism and its transformations. This model is made up of 8 differential equations. The following two formulae specify the behaviour of this model.

Irinotecan Exposure Control

In a cancer chronotherapy, an anticancer drug such as irinotecan is injected according to some control law over time. The control law can be represented by a series of parameterized events defining injection times and doses. An event is associated to the beginning and to the end of each injection. Parameters are used to characterize the lapse of time between consecutive injections. The injection control law is part of the system and represented as a component in its own right in the system.

The aim chosen here will be to minimize the toxicity (i.e., DNA damage on healthy cells, that are synchronized) while maintaining a fixed efficacy (the cancer cells lack circadian synchronization and thus efficacy will be supposed constant when the total amount of Irinotecan injected is constant, which is the assumption that the clinicians we collaborated with made). The parameters will thus be those defining a periodic step function with a fixed total area. Note that since the model considered only focusses on cellular pharmacokinetics (PK) and pharmacodynamics (PD), and not full-body PKPD, the law to be optimized is the exposure law instead of the injection law. Optimization of the injection law would follow a similar procedure but for a model with a defined target tissue and the corresponding PKPD.

Coupled Model Specification

Model Alignment

The first step of model coupling is model alignment for putting the models in the same format and normalizing molecule names. SBML versions of the irinotecan and p53/Mdm2 modules being available, they were imported in Biocham. The renaming of the variable representing DNA-damage was the only modification necessary in this precise case. More generally it would be necessary to rely on existing databases and ontologies to match corresponding entities in different models.

For the other models, we looked in parallel at the corresponding set of ordinary differential equations and at the available diagrammatic notation to write a set of Biocham reaction rules. Since ODEs can be automatically CPT11

DNA damage p53

Block cell cycle at checkpoints extracted back from the reactions, one can easily check that the reaction rule models are indeed coherent with the original ODE ones.

Apoptosis recover

Structural Coupling

The literature provides information about known structural links between the different building blocks to assemble them and compose the coupled model. Before that, let us examine the expected behaviour of the cell which is graphically depicted in Figure 5. Injections of irinotecan (CPT11) induce DNA damage. In reply to this, the cell reacts by activating protein p53, which blocks the cell cycle at a checkpoint. This arrest aims at repairing critical damage before DNA replications occurs, thereby avoiding the propagation of genetic lesions to progeny cells. Thus, while the cell cycle is arrested, the protein p53 will activate the DNA-damage repair mechanisms. If it is possible for the cell to recover, the cell cycle will be restarted; otherwise, if the damage is too extensive, the cell will undergo apoptosis.

As remarked in Section 3, literature provides evidence for the fact that, if a cell is exposed to irinotecan during the S phase of the cell cycle, then more DNA damage will be caused with respect to the other phases of the cell cycle [START_REF] Ohdo | Cell cycledependent chronotoxicity of irinotecan hydrochloride in mice[END_REF][START_REF] Zhou | Transcriptional regulation of mitotic genes by camptothecin-induced dna damage : Microarray analysis of doseand time-dependent effects[END_REF]. Keeping this fact in mind, we provided a characterization of the S phase in terms of the concentration level of CycA/Cdk2 (CycA for short) and we inserted in the irinotecan model a dependence from the S phase of the kinetic parameter involved in the production of DNA damage generated by the ternary reversible complexes SN38-Top1-DNA (Top1cc for short): such a parameter assumes a high value during DNA replication and a low value out of synthesis. In this way we linked the cell cycle model to the irinotecan one.

The structure of the coupling of the five models together is illustrated in Figure 6.

The link between the irinotecan model and the p53/Mdm2 one is given by DNA damage. In fact, irinotecan exposure causes DNA damage, which in turn triggers the activity of protein p53, that tries to recover DNA damage. The link between the cell cycle and circadian clock models comes from the experiments of [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF] and is reflected through a direct influence of CLOCK-BMAL1 (Bmal1) on the synthesis of Wee1, a kinase that delays or prevents entry into mitosis by phosphorylation of the Cdk1/CyclinB complex. This link uses the same structure as [9] since the Circadian clock model is the same. [START_REF] Zámborszky | Computational analysis of mammalian cell division gated by a circadian clock: Quantized cell cycles and cell size[END_REF] relied on a slightly different coupling that also modified, for unclear reasons, the reaction of CyclinB synthesis, whereas the aim here is to search for a coupling as simple as possible and satisfying the specification. Note that experimental results direct at a G2/M-transition focussed coupling but that for these experiments the cell-cycle model considered, even if it displays the four different phases, is centered around the restriction point following G1/S.

Bmal1 is also involved in the transcription of Top1 [START_REF] Yang | The molecular mechanism regulating the autonomous circadian expression of topoisomerase i in nih3t3 cells[END_REF]: this provides a link between the circadian clock and irinotecan models.

In order to link the p53/Mdm2 and cell cycle models, we inserted in the p53/Mdm2 model a rule which fixes that p53 activates p21, and two further rules imposing that p21 inhibits CycA and CycE, respectively. It is worth noting that we also investigated the possibility to abstract the previous expanded rules by letting p53 directly inhibit CycA and CycE. In the following, we will refer to this last version of the link as to the contracted one.

LTL(R) Specification of the Coupling

In this section, we show how the integration of temporal logic constraints and parameter optimization techniques can be used to compute kinetics for the coupled model. It is worth noticing that for this purpose, one can take advantage of LTL(R) formulae to express numerical constraints in a much more flexible way than by curve fitting, especially for oscillation constraints for instance.

The state transition structure is constituted of a simulation trace over a time window of 100 hours, containing the values of the system's variables of their first derivatives at discrete time points obtained by numerical integration (using Rosenbrock's implicit method for stiff systems).

We directly considered the model made up by the five components and all the linking rules, as illustrated in Figure 6, to perform the parameter research. For the sake of clarity, we will separately introduce each linking rule and the corresponding specification, but as a matter of fact we executed Biocham's parameter optimization procedure only ones to infer the unknown kinetic parameter values leading to the satisfaction of the conjunction of all the formulae given in the following.

The link between the circadian clock and irinotecan models (see Figure 6) has been encoded by means of the following reaction rule, that specifies a mass action law kinetics with parameter kbmaltop for the synthesis of Top1: MA(kbmaltop) for =[Bmal1 nucl]=>TOP1. The irinotecan model already included the following rule for the synthesis of Top1: top1 for =>TOP1. To keep the Top1 production limited and to constrain the correlation between the concentration values of Top1 and Bmal1, suitable values for top1 and kbmaltop such that property F1 holds have been found out. F1: Top1 is always lower than 1.5 and, whenever Bmal1 gets over 1 (before 85 time units), there exists a future state where Top1 is greater than 1.

LTL(R) : G([T OP 1] < 1.5 ∧ ([Bmal1 nucl] > 2.5) ∧ T ime < 85 → F([T OP 1] > 1)
). Results: we found out that the values top1=0.212 and kbmaltop=0.207 make F1 true.

The following Biocham rule encodes the link between the circadian clock and cell cycle models:

(ksweemp+ksweem*[Bmal1 nucl])/(kweem+kwpcn*[PER nucl-CRY nucl]) for =[Bmal1 nucl]=>Wee1.
While the cell cycle compounds oscillate with a period of approximately 23 hours, the circadian compounds exhibit a period close to 24 hours. To make the cell cycle properly be entrained by the circadian cycle, the cell cycle compounds have been required to oscillate with a period of approximately 24 hours, that is, we searched for values for the kinetic parameters involved in the above reaction rule so that property F2 is satisfied (see Figures 7 and8). circadian clock, is conserved when the p53/Mdm2 module is added but there is no irinotecan exposure. F3: Within a time interval of 100 time units, CycA is greater than 2.7 in at least 4 oscillations. LTL(R) : oscil([CycA], 4, 2.7). Results: suitable parameter values are k5321=0.487, kA21=0.00507, and kA53=0.283. Property F3 also turned out to be true when there is exposure to irinotecan but the p53/Mdm2 model is not taken into account. In fact, as expected, even if DNA damage occurs, when protein p53 does not act, the cell cycle is not affected, and thus CycA exhibits a regular oscillating behaviour.

Finally, to link the cell cycle and irinotecan models, the following rule has been considered MA(kdam) for TOP1cc=>DNAdam. already included in the irinotecan model and we made the parameter kdam depend from the S phase: it assumes a first value v1 during replication, a second value v2 out of replication, where v2 > v1. We searched for suitable values v1 and v2 for kdam so that the next property holds. F4: Whenever Top1cc gets above 0.2, there exists a future state when the first derivative of DNAdam gets above 0.15.

LTL(R) : G([T OP 1cc] > 0.2-> F(d([DN Adam])/dt > 0.15)).
Results: the property is verified for v1 equal to 1.42 and v2 equal to 1.89.

The combination of temporal logic constraints and parameter optimization techniques can also be used to validate the resulting model. As an example, the next property ascertains that, in case of repeated irinotecan exposure (and thus of sustained DNA-damage) the oscillations of CycA are affected. F5: When there is sustained DNA damage (after an initial period), the amplitude of CycA decreases before 73 time units and then stays low.

LTL(R) : F((T ime < 15) ∧ G([DN Adam] > 0.4)) → F((T ime < 73) ∧ G([CycA] < 2.

15)).

Results: with the expanded version of the links the amplitude of oscillations gradually decreases, satisfying the property. With the contracted one, oscillations are very irregular, as graphically depicted in Figure 9 (bottom panel).

Evaluation of the Parameter Search Procedure

The method used in Biocham to optimize parameter values with respect to LTL(R) properties consists in computing a continous satisfaction degree in [0, 1] for a temporal logic formula on a given simulation trace [START_REF] Rizk | On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology[END_REF], using an algorithm for computing validity domains of LTL(R) constraints with free variables instead of constants [START_REF] Fages | From model-checking to temporal logic constraint solving[END_REF][START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF] and then using the continuous satisfaction degree as fitness function for a continous optimization method.

Biocham uses the state-of-the-art nonlinear optimization method of Hansen and Ostermeier [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] named Covariance Matrix Adaptation Evolution Strategy (CMA-ES). A population of new candidate solutions is sampled according to a multivariate normal distribution of the parameters. The covariance matrix adaptation is a method to update the covariance matrix of this distribution. This method is a generalization of the approximate gradient and Hessian of a quasi-Newton method to an evolutionary algorithm for optimization problems with a black box fitness function on which no assumption is made. CMA-ES performs parameter search given an initial solution, stop and restart criteria, and a given search space. The search stops either when a given number of violation degrees have been computed or when the violation degree gets below a given threshold.

We searched for parameter values satisfying all F1 to F5 properties, each property being evaluated for a given set of models. The overall fitness of a parameter values set is the sum of the fitness of all these properties. This computation is done in parallel as well as the computation of the fitness of the population of solutions defined by CMA-ES. One 100h simulation of the complete model takes about 100s on a 3GHz processor and thus the evaluation of the 5 properties combined can take up to 500s. It took around 1000 evaluations of these five properties combined to find a satisfactory solution.

The execution time was 4 hours on 64 3GHz cores. Such a temporal logic constraint approach proved to be effective, allowing us to express relevant biological properties of the model (and concentration values that make specifications true) that could not be easily encoded as curve fitting problems for instance. This is the case in Figures 7 and8 which depict the behavior of the cell cycle when it is respectively entrained and not entrained by the circadian clock. While in the first case a period of approximately 23 hours is exhibited by CycA and CycB, in the second one the two compounds assume the same period of the circadian cycle, that is, approximately 24 hours. In Figure 7, the disruption of the first oscillation is due to the knock out of the entrainment reaction.

The set of the linking parameter values learned in Biocham with this procedure, together with the temporal logic formulae used for learning them, are recapitulated in Table 1

Optimal Control of Drug Exposure

The properties of this subsection deal with the control laws of irinotecan exposure. In order to deal with chronotherapeutics optimization for healthy cells while maximizing efficacy for tumor cells, we aimed at finding irinotecan exposure times and maximum amount that maintain toxicity low for healthy cells. More precisely, we modeled irinotecan exposure as rectangular boxes and we looked for maximum irinotecan quantity and for first exposure time, interval time between consecutive exposure and boxes width and height that keep DNA damage below a given threshold.

Evaluation of pulsatile exposure

In Figure 9 we show the behavior of the p53/Mdm2 DNA damage repair module when exposure is repeated every 24 hours. The plot puts in evidence how DNA damage increases after every exposure period. The oscillating trend of proteins p53 and Mdm2 is well highlighted. Furthermore, it is possible to notice the irregular behaviour assumed by CycA after exposure to irinotecan if the contracted link is used (bottom panel).

Optimization of the drug exposure law

To find the most efficient exposure law, we searched for the optimal schedule and the maximum amount of irinotecan such that DNA damage remains below a given threshold.

F6: DNAdam is always lower than 1 and total irinotecan exposure is greater than 50. LTL(R) :G([DN Adam] < 1) ∧ totalinjection > 50. We searched for parameter values that make F6 be satisfied with the lowest error, i.e., the values that maintain DNA damage below 1 and that minimize the distance between total irinotecan exposure and value 50. To avoid too short injections, the minimal injection length has been set to 1. Results: the maximum irinotecan exposure maintaining DNA damage low, are rectangular boxes with a width 1 (e.g., the lower bound we set to injection length) and a height of approximately 7. The first exposure should happen 23h30 after the initial state (chosen with CycA very low, i.e., in G1 phase), and a new exposure every cell cycle oscillation should then be done (see Figure 10, where the macro KDAM delineates the synthesis phase of the cell cycle).

It is worth noting that the formation of the ternary reversible complexes SN38-Top1-DNA (Top1cc) responsible for DNA damage follows the irinotecan exposure by a few hours and that the choice of irinotecan exposure described above corresponds at having Top1cc peaks out of the synthesis phase (remember that the production of DNA damage from Top1cc is higher during the S phase). On the other hand, the presence of Top1cc peaks during the replication phase leads at maximizing toxicity. Figure 11 shows the effect of the same exposure than for Figure 10 but with a 12h phase shift. With this phase shift, which can be attained for unsynchronized cells, DNA damage attains 1.7, that is a 70 percent increase compared to synchronized cells.

The next specification regards the DNA repairing power of the cell. F7: After an exposure to irinotecan is performed, DNA damage is able to go under the threshold of 0.1 before the next exposure.

LTL(R) : G(([CP T 11] > d) ∨ (([CP T 11] ≤ d)U([DN Adam] < 0.1))).
, where d depends on the dose of irinotecan. Before testing the property, we decided to parameterize the lapse of time between consecutive irinotecan exposures. Then we took advantage of the procedure learn parameters to find the minimum k such that, if one 10units-exposure is performed every k hours, then property F7 is true. Results: we found out that the minimum k multiple of 12 which makes F7 true is 36. Thus, one exposure every 36 hours should be performed in order to allow DNA damage to be recovered before the next exposure. Then we tried to see what it happens if, at each exposure, we double the irinotecan dose, that is, we expose to 20 units. In this case, one exposure every 48 hours should be done.

The last property requires the oscillating trend of proteins p53 and Mdm2 to stop before a new exposure. As for the previous specification, we parameterized the lapse of time between consecutive irinotecan exposures and we used the procedure learn parameters.

Results: the minimum k multiple of 12 which makes F8 true is 48.

Model Predictions for Circadian Clock Genes Knock-outs

Setup

Hereafter we describe how the cell cycle reacts to circadian gene/protein mutations in our coupled model.

This really amounts to verifying the predictive power of the model since, as already explained, the circadian entrainment is focussed on the G2/M transition, whereas the cell cycle model is focussed on the restriction point.

We explore what happens when a given compound is missing, that is, its concentration equals zero. To this aim, it is possible either to set the compound synthesis at zero, or to make the compound be absorbed by a "super-inhibitor" (e.g., the knock-out of a given compound C can be modeled by inserting in the model the rule Inhibitor + C → Inhibitor-C, where the initial concentration of Inhibitor is very high). As a matter of fact, both the alternatives have the same impact on the behavior of the coupled model. The mutations we take into consideration concern the mRNAs mPER, mCRY, and mBmal1. The simulations we provide in the following are up to 100 hours.

mPER=0. As shown in Figure 12, in this case the cell cycle period becomes bigger (approximately 28.5 hours), that is, the cellcycle is slowed down. Furhermore, the mean value of Wee1 is higher with respect to the one of the wild type phenotype. mCRY=0. In this case the behavior of the model is approximately the same of the previous one (see Figure 12). mBmal1=0. As illustrated in Figure 13 in this mutant the mean value of Wee1 is lower with respect to normal conditions and the cell cycle period is slightly smaller (approximately 23 hours).

Comparison with the Literature

In the following we itemize some facts we found in literature concerning the dependence of the cell division cycle on circadian rhythmicity/mutations and, when possible, we discuss the consistency with our results.

• The expression of several mammalian cell-cycle genes, including c-myc, Cyclin-D1, and mWee1, is regulated in a circadian manner [START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF]. As expected we also observe the circadian entrainment of cell-cycle genes in our in silico model.

• Overexpression of PER1 leads to apoptosis whereas inhibition of PER1 inhibits apoptosis. It appears that PER1 antagonizes the cell cycle in an oscillatory fashion similar to the manner in which it antagonizes the function of Clock-Bmal1 [START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF]. According to our experiments, a PER inhibition produces an increase of Wee1, an thus a mitosis inhibition. Note however that we observe a lengthened period and not a complete stop of the mitosis.

• In CRY deficient cells, the circadian rhythmicity is lost [START_REF] Gijsbertus | Mammalian cry1 and cry2 are essential for maintenance of circadian rhythms[END_REF], Wee1, over-expressed, and CyclinB, less active, loses rhythmicity [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF]. The effect on Wee1 and CyclinB is roughly consistent with our results.

There are also some KOs that were not directly comparable with our results since our model does not incorporate yet detailed DNA-damage pathways with ATR/ATM, Chk1/2 or cMyc:

• PER1 and TIM seem implicated on the DNA-damage response because both can be found complexed with the ATM and ATR kinases and the checkpoint kinases Chk2 and Chk1, respectively [START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF].

• The oscillatory expression of c-myc is abolished in mPER2 mutant mice, which could then result in an alteration of the p53 function [START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF].

Gene KOs Conclusion

We observe that for most of the knock-outs, the results of our coupled model are in accordance with experimental data, which considering the very simple specification used for the coupling is a quite interesting result.

For other mutations that should result in a complete stop of the mitosis, the result does not agree with the data since our model exhibits a slow down of the cell cycle but not an arrest in mitosis. This points out a weakness of the mammalian cell cycle model we have used. It is indeed driven by a constantly growing mass variable and focusses on the restriction point with few details on the G2/M transition. While this control of the mitosis by the mass variable is realistic in yeast, it limits the possibility of controlling the cell cycle in mammalian cells, it is thus virtually impossible to block at the corresponding checkpoint, even with a strong circadian coupling. These considerations motivate the use of cell cycle control models independent of the mass variable [START_REF] Gérard | Temporal self-organization of the cyclin/cdk network driving the mammalian cell cycle[END_REF][START_REF] Qu | Dynamics of the cell cycle: checkpoints, sizers, and timers[END_REF] allowing for more accurate predictions in this respect [START_REF] De | Model-based predictions of the influence of circadian clock genes knock-outs on the cell cycle[END_REF].

This evaluation of the model predictions on gene knock-outs also show that model-checking and parameter search are useful at the prediction stage: not finding any satisfactory parameter set when trying to strengthen the coupling in order to agree with the experimental result, indeed reveals a weakness in the structure of the individual models, which needs be revised in order to make the specification satisfiable.

Conclusion and Perspectives

In this paper, we have presented a coupled model of the mammalian cell cycle, circadian clock, p53-based DNA-damage repair, irinotecan intracellular PK/PD, and irinotecan exposure control, in order to study the influence of irinotecan drug in cancer chronotherapies. The coupling of the composite models has been achieved in Biocham using an original method based on LTL(R) temporal logic constraint solving, for representing the expected behavior of the coupled system, and on a continuous optimization evolutionary algorithm for inferring the values of the unknown coupling kinetic parameters of the models, as well as the exposure control parameters.

The maximization of antitumor effects and the minimization of the toxicity on healthy cells is the aim of any cancer therapy. The rationale of irinotecan chronotherapies is its toxicity on the cells in S phase only, the synchronization of the cell cycle by the circadian clock in healthy tissue cells, and the circadian disruption in mutated cancer cells. The resulting coupled model provides a valuable tool to investigate the drug influence on the cell cycle, reveal some weaknesses in the models, and ultimately infer some properties concerning the drug therapy and optimal exposure times and doses.

The predictive power of the coupled model was tested with respect to a limited set of mutants of the circadian clock genes. In the case of genes knock outs, we succeeded in considering temporal logic constraints over different traces corresponding to the mutations of different genes, that is, the initial condition of the trace relative to the knock out of a given set of genes is characterized by setting at 0 the parameters involved in the synthesis of the genes.

Although preliminary, the results obtained are very encouraging for our coupling method. In particular they showed that mass-entrained models of the cell-cycle have a limited possibility of entrainment by the circadian molecular clock. This motivates the use of non mass entrained cell cycle models like [START_REF] Gérard | Temporal self-organization of the cyclin/cdk network driving the mammalian cell cycle[END_REF][START_REF] Qu | Dynamics of the cell cycle: checkpoints, sizers, and timers[END_REF][START_REF] De | Model-based predictions of the influence of circadian clock genes knock-outs on the cell cycle[END_REF] which should not suffer from this limitation. The results also showed that the p53-Mdm2 DNA damage repair model of [START_REF] Ciliberto | Steady states and oscillations in the p53/mdm2 network[END_REF] should be improved in order to introduce a threshold above which the DNA is no longer repaired and the cell enters apoptosis. Last but not least, a PK/PD model of irinotecan in the body is missing to link the irinotecan injection law to the cell exposure model and optimize the drug injection law directly.

of Sciences of the United States of America, 100 [START_REF] Gérard | Temporal self-organization of the cyclin/cdk network driving the mammalian cell cycle[END_REF]

Figure 1 :

 1 Figure 1: Simulation plot of the cyclin concentrations during the mammalian cell cycle.

Figure 2 :

 2 Figure 2: Simulation plot of the mammalian circadian clock genes expression.

F 2

 2 p53 : Sustained DNA damage causes at least one oscillation of proteins p53 and Mdm2. LTL(R) :G(([DN Adam] > 0.2) → F(oscil([p53], 1)∧F(oscil([M dm2], 1)))).

F 3

 3 p53 : p53 oscillations are alternated by Mdm2 ones. LTL(R) : G(oscil([p53], 1) → X((¬oscil([p53], 1))U(oscil([M dm2 :: n], 1)))).

Figure 3 :

 3 Figure 3: Simulation plot of the P53/Mdm2 DNA-damage repair system.

F 1

 1 irin : In case there is no irinotecan, DNA damage equals 0. LTL(R) :G([CP T 11] = 0) → G([DN Adam] = 0). F 2 irin : If the irinotecan concentration is greater than 10, then there exist a future state when DNA damage exceeds the value 0.7 and then stays high. LTL(R) :G([CP T 11] > 10) → FG([DN Adam] > 3.5).

Figure 4 :

 4 Figure 4: Similuation plot of irinotecan metabolism.

Figure 5 :

 5 Figure 5: Schematic behaviour of the coupled model.

Figure 6 :

 6 Figure 6: Global schema of the coupled model.

Figure 7 :

 7 Figure 7: Simulation plot of the cell cycle (CycA, CycB) and circadian clock (mPER) with entrainment knock-out (no coupling). The cell cycle exhibits a free period of 23 hours.

F2:

  The period of CycA and CycB is 24. LTL(R) : period(CycA, 24) ∧ period(CycB, 24). Results: the values we found are ksweemp=0.521, ksweem=0.5, kweem=1, and kwpcn=2.Hereafter the Biocham rules introduced to link the p53/Mdm2 and cell cycle models are reported: MA(k5321) for =[p53]=>p21. MA(kA21) for CycA=[p21]=> . MA(kA21) for CycE=[p21]=> . As for the contracted version, the encoding is the following one: MA(kA53) for CycA=[p53]=> . MA(kA53) for CycE=[p53]=> . Again, suitable parameter values for k5321 and kA21 (kA53 in the second case) have been searched so that property F3, that expresses the CycA oscillating behaviour exhibited by the cell cycle model when entrained by the

Figure 8 :

 8 Figure 8: Simulation plot of the entrainment of the cell cycle by the circadian clock through coupling on Wee1. The period of the cell cycle is 24 hours.

Figure 9 :

 9 Figure 9: Simulation plot of DNA damage under pulsatile exposure to irinotecan every 24 hours with the p53/Mdm2 module. In the bottom panel, the contracted link is used, which results in very irregular CycA oscillations.

Figure 10 :

 10 Figure 10: Maximum exposure preserving DNA damage under threshold 1.

Figure 11 :

 11 Figure 11: DNA damage produced on phase-shifted cells with the same exposure law as in Figure 10.

F8:

  When exposed to irinotecan, p53 and Mdm2 are in a steady state, that is, their derivatives approach 0.LTL(R) : G(([CP T 11] > d) → ((d[p53] ≤ 0.05) ∧ (d[p53] ≥ -0.05) ∧ (d[M dm2 :: n] ≤ 0.05) ∧ (d[M dm2 :: n] ≥ -0.05)))., where d depends on the dose of irinotecan.

Figure 12 :

 12 Figure 12: Simulation of the coupled model with mPER=0.

Figure 13 :

 13 Figure 13: Simulation of the coupled model with mBmal1=0.

  (k3p+k3*[Cdc20])*(1-[Cdh1])/(J3+1-[Cdh1]) for _=[Cdc20]=>Cdh1. V4*[Cdh1]/(J4+[Cdh1]) for Cdh1=>_. epsilon*k11p for _=>Cdc20_T. MA(epsilon*k11) for _=[CycB]=>Cdc20_T. MA(k12) for Cdc20_T=>_. k13*[IEP]*([Cdc20_T]-[Cdc20])/(J13+[Cdc20_T]-[Cdc20]) for _=[IEP]=>Cdc20. (k14/(J14+[Cdc20])+k12)*[Cdc20] for Cdc20=>_. epsilon*k33 for _=>PPX. MA(k34) for PPX=>_. k31*[CycB]*(1-[IEP])/(J31+1-[IEP]) for _=[CycB]=>IEP. k32*[PPX]*[IEP]/(J32+[IEP]) for IEP=[PPX]=>_. k27*mass*(if Rb_hypo/Rb_T >0.8 then 0 else 1) for _=[massT]=>GMT. MA(k28) for GMT=>_. epsilon*mu*nbcells*[GMT]for _=[GMT]=>massT.

Table 1 :

 1 Parameter values learned in Biocham, values found, and temporal logic formulae used for learning them.

	Parameter Value Formula
	kbmaltop	0.207	F1
	top1	0.212	F1
	ksweemp	0.521	F2
	ksweem	0.5	F2
	kweem	1.12	F2
	kwpcn	5	F2
	k5321	0.486	F3
	kA21	0.00507	F3
	kA53	0.283	F3
	v1	1.42	F4
	v2	1.89	F4

  :14806-14811, December 2003.

	MA(V6)	for _=[CycE-Kip1]=>CycE.
	(MA(k25),MA(k25r))	for CycE+Kip1<=>CycE-Kip1.
	MA(V6+V8)	for CycE-Kip1=>_.
	epsilon*k29*E2F_A*mass	for _=[massT]=>CycA.
	MA(k30)	for CycA=[Cdc20]=>_.
	MA(V6)	for _=[CycA-Kip1]=>CycA.
	(MA(k25),MA(k25r))	for CycA+Kip1<=>CycA-Kip1.
	MA(V6)	for CycA-Kip1=>_.
	MA(k30)	for CycA-Kip1=[Cdc20]=>_.
	epsilon*k5	for _=>Kip1.
	MA(V6)	for Kip1=>_.
	MA(k10)	for _=[CycD-Kip1]=>Kip1.
	MA(V8)	for _=[CycE-Kip1]=>Kip1.
	MA(k30)	for Cdc20+CycA-Kip1=>Kip1+Cdc20+CycA-Kip1.
	k22*E2F_T	for _=>E2F.
	MA(k22+k23p)	for E2F=>_.
	MA(k23)	for E2F=[CycA]=>_.
	MA(k23)	for E2F=[CycB]=>_.
	epsilon*k1p	for _=>CycB.
	epsilon*k1*([CycB]/J1)^2/(1+([CycB]/J1)^2)	for _=[CycB]=>CycB.
	MA(V2)	for CycB=>_.

http://contraintes.inria.fr/supplementary material/TCS-CMSB09/.

Acknowledgements

This work was supported by the EU FP6 STREP project TEMPO on cancer chronotherapies and is now supported by the ERASysBio project C5Sys concerning circadian and cell cycle clock systems in cancer. We acknowledge fruitful discussions with the partners of this project, in particular with Francis Lévi, Jean Clairambault and Annabelle Ballesta.