Institut Camille Jordan

I would like to thank my advisors Mr. Vitaly Volpert and Mr. Samuel Bernard for the time they dedicated to my work. I thank them for their trust and for letting me develop my own ideas, helping me to pave my way in a theme that is still open for many questions. Their remarks and advice were very helpful and enriching. I am grateful for all the knowledge and modeling techniques they transmitted to me, our discussions were always bringing new insights and clarifying confusing things. I did not feel along the thesis period as a student obliged to do his work, Mr. Volpert and Mr. Bernard gave me the trust to conduct with responsibility my own project; I appreciate this too much! I would like to thank Mr. Marc LeFranc and Mr. Didier Gonze for the interest they showed to my work and for accepting to review my thesis. I thank also Mme Angélique Stéphanou, Mr. Daniel Le Roux and Mr. Dominique Barbolosi for accepting to be part of my thesis examination jury. I thank a lot Mme Elizabeth Mironescu the head of the Institut Camille Jordan and the administrative staff for all their help and support.

A lot of thanks for the Dracula team at INRIA. I thank first Mr. Mostafa Adimy the head of the team for accepting me as a member of Dracula. I thank Thomas for all the discussions we had, and especially for the leukemia project that he let me participate to. I thank Laurent for his support during the ATER campaign. I am grateful a lot for Fabien for all the time and discussions we had, his support and advice were so important and encouraging. I thank finally the team secretary Mme Caroline Lothe for all the help with the administrative procedures. I thank my cousin Nader for the lot of support he and his wife gave me when I first arrive in France. I thank him also for the great time I spent at LAU working on the Ceder Project. I Thank a lot Doron Levy for inviting me to CSCAMM at the University of Maryland to work on CML modeling. It was a great opportunity working with him. I am Grateful for that. I thank also a lot Geoffry Clapp with whom I worked on CML. Getting into the work with him was very fast and efficient. I thank him a lot for the pleasant time and discussions we had. I thank all the PhD students that I met during my thesis. Special thoughts go for my colleagues at the Insitut Camille Jordan, namely Ivan, Evrad, Adrianne, Hassan; Also my colleagues at INRIA, Abdelnasser, Loic and Marine. Many thanks also go to my Lebanese friends especially Abbas for all the time spent together and for the support. This work could not be done without the help and support of my wife whom I cannot thank enough. She is the person who was always there for me in the bad times as in the good times, listening to my complaints, easing the lot of stress and frustration I had during my thesis.

Finally I thank my parents and my brother for their infinite support; they are the reason why I had accomplished this work. They are a God's blessing. I dedicate this thesis to them.

List of publications

Published papers

• R. El Cheikh, S. Bernard and N. El Khatib. Modeling circadian clock-cell cycle interaction effects on cell population growth rates (2014), J. Theor. Biol., 363:318-331.

• R. El Cheikh, T. Lepoutre

Submitted papers

• Geoffrey D. Clapp, Thomas Lepoutre, Raouf El Cheikh, Samuel Bernard, Jérémy Ruby, Franck E. Nicolini, Doron Levy, BCR-ABL transcripts variations in chronic phase chronic myelogenous leukemia patients on imatinib: Possible role of the autologous immune system.

Papers in preparation

• R. El Cheikh and S. Bernard: Multi-scale model for the regulation of the cell cycle by the circadian clock.

• S. Bernard, R. El Cheikh, H. Huttner, O. Bergmann and J. Frisèn. A mathematical model to estimate the age of meningioma.

Abstract

This thesis is dedicated to the development of a multiscale mathematical model that describes the regulation of the cell cycle by the circadian clock. What motivated this work is the fact that several tumorigenic diseases are linked to circadian rhythms disruption. We would like to understand the effect of circadian rhythms on the proliferation of a cell population and hence give plausible explanation for diseases that arise from circadian clock disruption.

The mammalian cell cycle and the circadian clock are two molecular processes that operate in a rhythmic manner and exquisite precision. On one hand, the cell cycle is driven by the rhythmic activity of cyclin-dependent kinases which dictate the time a cell must engage mitosis and the time it must divide giving birth to two daughter cells. On the other hand, the circadian clock is a system of transcriptional and translational feedback-loops that generates sustained oscillations of different mRNAs and proteins with a period of approximately 24 h. It turns out that several components of the circadian clock regulates various cyclin-dependent kinases at different stages of the cell cycle. This makes the circadian clock a key player of the temporal organization of the cell cycle and makes these two biological processes act as two tightly coupled oscillators.

Our modeling approach consists of using a molecular-structured partial differential equation that describes the proliferation of a cell population. Proliferation depends on the coupled cell cycle-circadian clock molecular state of cells. Due to the large number of molecular components involved in the cell cycle-circadian clock system, the problem becomes of highdimensionality and specific numerical techniques are needed to solve the equation.

As a first step, we simplify the problem, and use a system of transport partial differential equations structured by the time spent by cells in a phase of the cell cycle. This system is coupled to the molecular one, via transition coefficients that depend on the molecular state of cells. Even tough it is a simplified version, this model has the novelty of combining both population and intracellular levels. We use it to study the entrainment of the cell cycle by the circadian clock and the effects of regulation on the net growth of cells.

Afterwards, we pass to the fully multi-scale model and use the particle method to circumvent the high-dimensionality aspect. This method consists of representing the population of cells by a large number of particles, each having its own set of properties, position and weight. These properties evolve in time according to a system of ordinary differential equations, so that the particles simulate individual cells evolving in the molecular state space.

The solution of the system can be reconstructed from individual particles and shown to converge to the solution of the initial system. The main advantage of particle methods is their usefulness for high dimensional problems where classical numerical methods as finite difference/volumes/elements fail.

The thesis contains three main chapters, a conclusion and an Appendix:

• Chapter 1

This chapter covers biological fundamentals (enough though for our mathematical study) of the circadian clock and the cell cycle and mathematical models used to study them. It contains three sections; in the first one, we give a brief review of the history of circadian concepts and the molecular mechanism of mammalian circadian clock. We then review the evolution of mathematical models used for the circadian clock, starting from phase response curve models to detailed molecular ones. The second section of this chapter is a description of the temporal organization of the cell cycle and its mathematical modeling. In the third section we introduce renewal and structured partial differential equations used to model cell proliferation. We review main theorems on the net growth of a proliferating cell population under circadian control. Part of this chapter, namely the presentation of the mathematical models for the circadian clock, and sections 2,3 are based on the paper • Chapter 2

In this chapter, which is based on our published paper R. El Cheikh, S. Bernard and N. El Khatib, Modeling circadian clock-cell cycle interaction effects on cell population growth rates (2014) J Theor Biol, 363:318-331.,

we introduce the mathematical model that represents the first step of our study; modeling the regulation of the cell cycle by the circadian clock. The novelty of this model is that it contains both molecular and population levels. We investigate the way the cell cycle entrains to the circadian clock with different rational period ratios and characterize multiple domains of entrainment. We show that circadian clock increases the growth rate in cells with autonomous periods of the cell cycle around 24 h and above 48 h. We study the effect of mutation of circadian genes on the growth rate of cells and show that disruption of the circadian clock can lead to abnormal proliferation. We obtain non-intuitive results on growth rates that could not be obtained by a population or a molecular model alone. Our model offers new insights on the influence of the circadian clock on the growth of a cell population.

• Chapter 3

In this chapter, we present a multiscale model for the regulation of the cell cycle by the circadian clock. It consists of a non-linear transport equation of the form ∂ t ρ(x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = L(x, λ)ρ(x, t, λ).

This equation describes the evolution of cells density ρ with a convective term u that depends on the molecular circadian-cell cycle system components x = (x 1 , • • • , x d ) and a growth term L that describes the population dynamics. The variable λ is a cellular state parameter that induces heterogeneity among cells, and ψ is a function that computes population statistical quantities such as average molecular concentrations. This makes the cell density molecularly structured and hence confers to it a high-dimensional aspect. One of the aims of this chapter is to adapt a numerical method to solve the transport equation and to develop a code for it (in C, given in the appendices) that can be used with different convective velocities u. The use of "structured" transport equations in biology is usually limited to few dimensions, hence our code can pave the way for new applications with high-dimensional equations. This chapter is divided into three sections; in the first one, we present the model structure, the additional modeling assumptions and details about the equations. In the second one, we introduce the particle method and explain its theoretical background and the way it is used to solve high dimensional equations. In the last part, we present a numerical test case to illustrate the particle method and to test our code and then present results obtained with our model; namely, results that confers to this model its multi-scale nature, like the heterogeneity among cells and its influence on the growth rate, the dependence of growth on the total cell number and the connectivity between cells and its implication on population synchronization.

• Conclusion

We finish the main part of the manuscript with a conclusion in which we summarize the work and results obtained during this thesis. We discuss possible improvements of our modeling approach and the continuation of this study.

• Appendix

The appendix contains mainly two sections. In the first on we present a study (made during the thesis period) that is not directly related to the topic of circadian clock-cell cycle coupling; however it still belongs to theme of age-structured equations and tumor modeling. We present a model that we developed to estimate the time of appearance for meningioma tumor for several patients. We had clinical data about the average age of tumor cells, their Ki-67 index and the tumor volume. Using minimization techniques, we fitted the data and gave an estimate for the tumor age. In the second section of the appendix, we give a commented version for the code source files that we developed about the particle method. In addition, the appendix contains some supplementary materials as table values and additional formulas that we used in this manuscript.

Résumé

Cette thèse est dédiée au développement d'un modèle mathématique multi-échelle pour la régulation du cycle cellulaire par l'horloge circadienne. Ceci est motivé par le fait que plusieurs études ont montré un lien direct entre certains cancers et un dysfonctionnement du mécanisme de l'horloge circadienne. Le but est de comprendre l'effet des rythmes circadiens et leur perturbation sur la prolifération d'une population de cellules.

Le cycle cellulaire et l'horloge circadienne sont deux processus moléculaires qui fonctionnent de manière rythmique et avec une précision exquise. D'une part, le cycle cellulaire est contrôlé par l'activité rythmique des cyclines CDK (cyclin-dependent kinases) qui dictent le temps pour lequel une cellule doit s'engager dans une mitose et le temps de division pour lequel une cellule donne naissance à deux cellules filles. D'autre part, l'horloge circadienne est composée d'un système de boucles rétroactives transcriptionnelles et translationnelles qui génèrent un régime d'oscillations permanent d'une multitude de messagers ARN et protéines avec une période de 24 h. Il s'avère que différents composants de l'horloge circadienne régulent plusieurs cyclines CDK lors de différentes étapes du cycle cellulaire. Ceci donne à l'horloge circadienne un rôle principal dans l'organisation temporelle du cycle cellulaire et fait de ces deux processus biologiques deux oscillateurs intimement couplés.

Notre approche pour décrire la prolifération d'une population de cellules consiste à considérer une équation de transport structurée par les contenus moléculaires du système couplé horloge circadienne/cycle cellulaire. En raison du nombre élevé des composants moléculaires intervenant dans ce système, le problème devient hautement multidimensionnel. Des méthodes numériques spécifiques sont requises alors pour la résolution. Dans un premier temps, pour simplifier le modèle, nous considérons un système d'équations de transport structurées par le temps passé par les cellules dans une phase du cycle cellulaire. Ce système est couplé au système intracellulaire via des coefficients de transition dans l'équation du transport qui dépendent du mécanisme moléculaire. Malgré le fait que ce premier modèle représente une version simplifiée, il a la nouveauté de décrire les deux niveaux; populationel et intracellulaire. Nous l'utilisons pour obtenir des résultats sur l'entraînement du cycle cellulaire par l'horloge circadienne, sur l'effet de la régulation et la perturbation de l'horloge sur la croissance cellulaire. Ensuite nous passons au modèle multi-échelle et nous contournons la difficulté de la haute-dimension de l'équation par l'utilisation de la méthode de particules. Cette méthode consiste à représenter la densité des cellules par un grand nombre de particules, dont chacune possède des propriétés héritées du système intracellulaire comme le contenu intracellulaire et les valeurs de paramètre. Ces propriétés évoluent selon un système d'équations différentielles ordinaires, de sorte que les particules simulent des cellules qui se déplacent dans l'espace de contenus moléculaires. La solution de l'équation peut être reconstruite de la distribution des particules et converge vers la solution classique. L'avantage principal des méthodes de particules est leur utilité pour la résolution des équations de haute dimension où les méthodes numériques classiques, comme différences/volumes finis, sont pratiquement inutilisables.

Cette thèse contient trois chapitres. Dans le premier, nous expliquons d'une façon simplifiée-suffisante cependant pour ce travail-la théorie biologique du cycle cellulaire et de l'horloge circadienne, et on passe en revue les différents modèles mathématiques utilisés pour les décrire. Ce chapitre est divisé en trois parties; la première est consacrée à l'horloge circadienne. Dans cette partie, nous racontons brièvement l'histoire des rythmes circadiens. Ensuite, nous expliquons le mécanisme moléculaire des ces rythmes chez les mammaires et nous retraons l'évolution de leurs modèles mathématiques, à commencer par les modèles de réponses de phase PRC (phase response curve) jusqu'au modèles moléculaires détaillés. Dans la deuxième partie nous décrivons l'organisation temporelle du cycle cellulaire et sa modélisation mathématique. Dans la dernière partie, nous introduisons les équations aux dérivées partielles de renouvellement qui sont utilisées pour la modélisation de la prolifération d'une population de cellules. Nous exposons les différents théorèmes liés à l'impact d'un contrôle circadien sur le coefficient de croissance.

Dans le deuxième chapitre nous introduisons le modèle simplifié qui étudie la régulation du cycle cellulaire par l'horloge circadienne. Ce modèle a la nouveauté de décrire la prolifération d'une population de cellules tout en tenant compte du système couplé horloge circadienne-cycle cellulaire. Dans ce chapitre, nous étudions comment le cycle cellulaire peut être entrainer par l'horloge circadienne sous multiples périodes rationnelles et nous caractérisons différents domaines d'entraînement. Nous montrons que le couplage augmente le taux de croissance pour des cycles cellulaires d'une période proches de 24 h et plus grande que 48 h. Ensuite, nous étudions l'effet des mutations des gènes de l'horloge circadienne sur le taux de croissance. Nous démontrons qu'une perturbation du mécanisme circadien résulte en une croissance anormale des cellules. Nous obtenons des résultats contre-intuitifs qui n'auraient pas pu être obtenus avec un modèle moléculaire ou populationnel seulement. Ceci montre que notre modèle offre des nouvelles perspectives de l'influence de l'horloge circadienne sur la prolifération cellulaire.

Dans le troisième chapitre nous introduisons le modèle multi-échelle. Ce modèle est constitué d'une équation non-linéaire de transport de la forme ∂ t ρ(x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = L(x, λ)ρ(x, t, λ).

Cette équation décrit l'évolution de la densité des cellules ρ avec un terme de convection u qui dépend des composants moléculaires x = (x 1 , • • • , x d ) de l'oscillateur couplé cycle cellulaire/horloge circadienne. La variable λ est un paramètre d'état qui induit une hétérogénéité parmi les cellules, et ψ est une fonction qui calcule des quantités statistiques au niveau populationel, par exemple, une concentration moyenne d'un composant moléculaire. Ceci rend l'équation de transport structurée par les contenus moléculaires et lui donne un aspect hautement multidimensionnel. Un objectif principal de ce chapitre est d'adapter une méthode numérique pour la résolution de cette équation, puis développer un code (en langage C) qui pourra être utilisé avec différentes vitesses de convection u. L'utilisation des équations de transport "structurées" dans des applications biologiques est souvent limitée à des dimensions non élevées. Ce code permettra d'élargir les domaines d'application à des équations de dimensions élevées. Ce chapitre est divisée en trois parties; dans la première nous introduisons le modèle en expliquant les différentes hypothèses de modélisation que nous y incluons. Dans la deuxième, nous présentons d'une manière simplifiée la théorie derrière la méthode de particules. Dans la dernière partie, nous résolvons un cas test numérique pour illustrer la méthode de particules et pour tester notre code. Ensuite, nous présentons les résultats du modèle, surtout ceux liés à sa nature multi-échelle, comme l'hétérogénéité de la population, la dépendance de la prolifération sur le nombre total des cellules et le paramètre λ, et l'implication de la connectivité intercellulaire sur la synchronisation des rythmes.

List of Figures

Chapter 1

State of the art 1.1 Introduction

The initial objective of this thesis was to study the impact of circadian rhythms on the dynamics of a population of proliferating cells, using multi-scale modeling tools. We would like to give an explanation for the disorders that may arise when these rhythms are disrupted. For that, three pieces should be put together to complete the whole picture: the circadian clock, the cell cycle and the proliferation of a cell population. These three pieces can be put together by understanding the way the cell cycle and the circadian clock are connected. Then by understanding how to link this coupling to the proliferation of a cell population. This chapter serves as the biological and mathematical background for our study. We give biological details about the cell cycle and the circadian clock and passes in review mathematical tools used to model them. This helps knowing what has been done until now in related research fields and what additional studies our work brings.

The first section focuses on the circadian clock, the daily rhythm that regulates physiology and behavior in our body. We present a short history of circadian clock research, starting from early botanic experiments to recent genetic discoveries. We then present the genetic mechanisms for the mammalian circadian clock, and describe the main transcriptionaltranslational feedback loops that generate sustained oscillations. We end this section with a review of mathematical models constructed to study the circadian clock. We look at phase response curve models, namely the works of Colin S. Pittendrigh and Jürgen Aschoff, and how they were used to study entrainment properties. We then introduce generic models, especially Goodwin-type models and explicit time-delay models. We discuss their properties and the importance of time delays for generating oscillations. We finally present molecular models, starting with basic models of this type that describe a single feedback-loop and finishing with recent complex models including detailed molecular processes involved in the circadian mechanism.

The second section is devoted to the cell cycle, its temporal organization, molecular mechanisms and its mathematical modeling. We describe how the cycle is driven by cyclindependent kinases and how a cell enters into mitosis and divides giving birth to two daughter cells. We explain the different stages a cell must pass through to accomplish a healthy cycle, portraying it as a succession of transitions that should occur in the right time and conditions. This yields a dynamical view of the cell cycle that helped mathematicians modeling it using differential equations. In the end of this section, we present some of these works, mainly the models developed by Tyson and Novak.

In the third section we introduce cell renewal equations and structured cell division models. We explain how these equations are used to model cell proliferation. We review the main results obtained for these models on the growth of a population under periodic forcing. This represent a crucial point for our work, since these results are among the first ones that deal with the implication of circadian control on cell proliferation.

Circadian clock 1.2.1 History of circadian rhythms

This subsection about the History of circadian rhythms is based on chapter 1, "A History of Chronobiological Concepts" in U Albrecht (ed.), The Circadian Clock from Protein Reviews edition volume 12, written by S. Daan [START_REF] Daan | A History of Chronobiological Concepts[END_REF].

Humans have always been aware of the ubiquity of rhythms governing their life. The alternation of seasons, day and night, the migration of birds, the periods of sowing and harvesting, plants leaves movements, wake and sleep phases, are all examples that drew the attention of natural philosophers and let them wonder if there might be a "force" orchestrating these rhythmic phenomena.

Botanists were especially curious about the daily movements of plants leaves. In the beginning, it was not clear that most of these rhythms came from within plants, and not from factors related to earth rotation, given that these movements are in great harmony with the alternation of everyday light and dark cycles. The first breakthrough was the experiment done by the French astronomer Jean Jacques d'Ortous de Mairan in 1729 who hid a mimosa plant in a dark cupboard and observed that the movement of its leaves persisted in the absence of light [START_REF] Mairan | Observation botanique[END_REF]. However, he was cautious about his conclusions and suggested that the plant still sense the sun without seeing it, and that other factors, such as temperature and humidity, might be still involved in leaves movements. In 1758 the French physician and botanist Henri-Louis Duhamel du Monceau reported leaf movements in a dark cave with constant humidity and temperature [START_REF] Duhamel | La physique des arbres[END_REF]. Later on, by the year 1875; the Swiss botanist August de Candolle added a great contribution to de Mairan and Duhamel experiments. He realized that movements of mimosa leaves, kept in constant light, followed a 2 h shorter cycle. This free running cycle is sufficient to claim the endogenous nature of rhythmicity. It would be unreasonable to attribute this cycle to any external factor since movements were no longer in synchrony with the surrounding environment.

Despite these convincing experiments, botanists were reluctant to conclude about the true endogenous nature of these rhythmicities. One of the skeptics was the famous German botanist Wilhelm Pfeffer. In his early career, he believed that there were other factors related to earth rotation that produced these rhythms. He argued that leaf movements tend to damp out in constant environment, which make the resulting free-running rhythm like a "Nachschwingungen" (after oscillations) phenomena [START_REF] Pfeffer | Die Periodischen Bewegungen der Blattorgane[END_REF]. As a consequence, for nearly 200 years and until the beginning of the twentieth century, botanists kept looking for hidden factors related to earth rotation that could explain the nature of these rhythmic processes.

During the twenties and early thirties, two botanists ended up this controversy and undoubtedly demonstrated the endogenous nature of plant movements. The Dutch botanist Antonia Kleinhoonte observed that after a shift in the light-dark cycle, the phase of the rhythm of Canavalia ensiformis in constant conditions had no longer any relationship to the day and night outside [START_REF] Kleinhoonte | Über die durch das Licht regulierten autonomen Bewegungen der Canavaliablätter[END_REF]. Following her work, the German biologist Erwin Bünning confirmed that plants, held in constant light, showed sustained rhythmic movements with a cycle length deviated from 24 h [START_REF] Bünning | Über die tagesperiodischen Bewegungen der Primarblatter von Phaseolus multiflorus. II. Die Bewegungen bei Thermokonstanz[END_REF]. Both drew the conclusion that this could not be attributed to any factor related to earth rotation and that it should be generated from within the plant.

Meanwhile, premises for the presence of internal rhythms in animals were unraveled. In 1900, the English zoologist Frederick William Gamble reported the first animal rhythm yielding evidence of a deviation from 24 h in pigmentation change in the crustacean Hippolyte varians [START_REF] Gamble | Hippolyte varians: A study in colour change[END_REF]. Some years later, the American physiologist, Curt Richter added to Gambles findings and observed the persistence of rhythmicity in rats activity in constant environmental conditions with a period shorter than 24 hours [START_REF] Richter | A behavioristic study of the activity of the rat[END_REF].

The proof of the endogenous nature of these rhythmic phenomena was a turnaround in this field and paved the way for new discoveries, especially during the fifties and sixties, when several new properties about biological rhythms were identified.

Scientists wondered whether these endogenous rhythms are innate or if they are, due to decades of evolution, a result of the eternal alternation of light and day cycles. To answer this query, Jürgen Aschoff, a German medical doctor, raised mice for several generations in constant light, and found that the daily rhythms of the last generation were as robust as those seen earlier, despite the persisting absence of a lightning cycle [START_REF] Aschoff | Tagesperiodik bei Mäusestämmen unter konstanten Umgebungsbedingungen[END_REF]. He also raised chickens (developing inside the egg under constant conditions) and found that they expressed a normal daily rhythm, indicating that presence of a light dark cycle was not needed for the development of rhythmicity [START_REF] Aschoff | Angeborene 24-Stunden-Periodik beim Käcken[END_REF]. These results were confirmed later on by Vijay Sharma in Bangalore, who grew Drosophila in constant light for more than 600 generations, and still found no loss of circadian rhythmicity [START_REF] Sheeba | Persistence of eclosion rhythm in drosophila melanogaster after 600 generations in an aperiodic environment[END_REF]. This led to the concept of innateness; which means that rhythmicity was not only endogenous, but also transmitted through the genetic material. This implies that the rhythmic environment is not essential for the function of these rhythms but plays an adjustment role.

All these results formed a solid foundation for the evidence of internal rhythmicity in living species and the word circadian saw the light and was first used by Franz Halberg to insist about the fact that "about a day" is the nature of these rhythms rather than 24 h [START_REF] Halberg | Physiologic 24-hour periodicity in human beings and mice, the lighting regimen and daily routine. Photoperiodism and related phenomena in plants and animals[END_REF]. The conclusion that most living species have developed an internal endogenous oscillating system was finally admitted and nowadays a free-running rhythm, with a cycle deviating from 24 h in constant conditions is considered as a sufficient proof of the endogenous generation.

After early experiments done by Antonia Kleinhoonte, advances in the biological theory of circadian rhythms promoted the emergence of multiple physical experiments and mathematical models to study further physical properties of circadian oscillations. Kleinhoonte showed that single stimuli of light can shift the phase a circadian rhythm and that a stimuli as brief as 1 min of light can shift the rhythm of the large leaves of the bean plant Canavalia ensiformis [START_REF] Kleinhoonte | Über die durch das Licht regulierten autonomen Bewegungen der Canavaliablätter[END_REF]. Following Kleinhoonte work, two persons, Colin S. Pittendrigh and Jürgen Aschoff, laid the foundation of a new field of study for circadian rhythms: entrainment; a concept that describes the way oscillators are synchronized by external signals. Aschoff introduced the concept of "Zeitgeber" (literally "time giver"), to indicate synchronization of a self sustained oscillator by an external signal and Pittendrigh introduced the concept of phase response curve PRC to describe the way a perturbation affect the oscillations period [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF]. They had different views on the way entrainment should be described and developed two theories: parametric and non-parametric entrainment. Aschoff and Pittendrigh collaborators, namely Rütger Wever and Arthur Winfree, widened the research in this field; Wever developed what is, perhaps, the first mathematical model of the circadian clock, to study the influence of light fluctuation on the circadian period [START_REF] Wever | Ein mathematisches Modell für die circadiane Periodik[END_REF]. Winfree did a marvelous job describing the way multiple oscillators synchronize together and characterized type 1 and type 0 response curves [START_REF] Winfree | Biological rhythms and the behavior of populations of coupled oscillators[END_REF][START_REF] Winfree | Integrated view of resetting a circadian clock[END_REF]. Based on these mathematical frameworks, a broad class of generic mathematical models were developed to study different properties of circadian oscillations.

Since most of early studies of circadian rhythms were carried out on plants, which do not possess a nervous system, the concept of circadian rhythms was perceived as a diffuse capacity of the entire organism. Surprisingly, several experiments discarded this idea. In 1960, Curt Richter did a series of brain lesions in rats and suggested that the hypothalamic area was a candidate for harboring a circadian pacemaker [START_REF] Richter | Sleep and activity: Their relation to the 24-hour clock[END_REF]. In 1972, Friedrich Stephan and Moore R.Y. observed that ablation of a region in the hypothalamus called the suprachiasmatic nucleus (SCN) caused rhythmicity to disappear. And by the year 1979; Inouye and Kawamura demonstrated that the SCN could be isolated from the surrounding brain tissue and then retains its rhythmicity in multiple unit activity, while the rest of the brain and the animal would become arrhythmic. From these experiments emerged the idea of a central pacemaker that orchestrates rhythms in the whole organism. The definite proof was given by Martin Ralph, who exploited the first mammalian circadian mutant and transplanted SCNs from homozygous mutant embryos into SCN-lesioned wildtype hosts and vice versa. He showed in this elegant experiment that the donor always determined the period of the restored circadian rhythm [START_REF] Ralph | Transplanted suprachiasmatic nucleus determines circadian period[END_REF].

The idea of pacemaker and entrainment properties motivated scientists to search for the input pathways by which the Zeitgebers (the light-dark cycle) were perceived. Junko Nishiitsutsuji-Uwo first demonstrated, by optic tract sections on the cockroach, that the compound eyes mediate the light information to the pacemaker [START_REF] Nishiitsutsuji-Uwo | Central nervous system control of circadian rhythmicity in the cockroach[END_REF]. Recently several experiments showed that cockroaches possess multiple photoreceptive inputs into their complex circadian systems, including the eyes, the pineal and other light-sensitive elements. This is not the case in mammals. Circadian entrainment by light employs an exclusive retina-hypothalamic pathway. Curt Richter demonstrated that blinding rats and monkeys abolishes entrainment by light while leaving their rhythm intact. Then, in the beginning of 1980s Groos and Mason showed that the visual field in the retina is connected to single cells in the SCN [START_REF] Groos | The visual properties of rat and cat suprachiasmatic neurones[END_REF]. This was confirmed two decades later by Russel Foster and his group; who found a specific network of intrinsic photosensitive ganglia in the inner retina. These ganglia contain pigment melanopsin and communicate information directly to the SCN [START_REF] Provencio | A novel human opsin in the inner retina[END_REF].

The fact that the circadian mechanism is encoded in the genetic material could not be made clear until the first discoveries about the DNA structure and genetic code. Ronald Konopka studied the aberrations in the circadian system controlling daily pupal eclosion in Drosophila. He identified three mutants: one with long period of 28 h (Per L), one with 20 h (Per S) and one arrhythmic (Per 0) [START_REF] Konopka | Clock mutants of drosophila melanogaster[END_REF]. The crucial step in this direction was the finding by a research group at Brandeis University that a protein, labeled PER, was rhythmically produced in the fly brain with a peak in the early night, while per RNA showed a similar expression pattern, about 6 h in advance of the protein. This led Paul Hardin, Jeff Hall and Michael Rosbash in 1990 to propose the idea of the transcription-translation feedback loop [START_REF] Hardin | Feedback of the drosophila period gene product on circadian cycling of its messenger RNA levels[END_REF]. The proposition was that the Per gene is transcribed in the nucleus of certain cells, leading to messenger RNA leaving the nucleus into the cytoplasm, and to translation into PER protein at the ribosomes, followed by return of the protein into the nucleus to suppress further transcription of the gene.

In mammals, the first mutation involved in the clock mechanism was detected by accident in Syrian hamsters, Mesocricetus auratus. The circadian period in DD in hamsters varies between 23 and 25 h. In 1986 Martin Ralph, depicted the remarkable actogram of one hamster displaying a circadian rhythm of activity with a period of 22 h. His insight to breed this individual led to the finding that the mutation followed simple Mendelian inheritance [START_REF] Ralph | A mutation in the circadian system in golden hamsters[END_REF]: A single gene mutation (called tau) must be involved. Homozygous mutants had a period around 20 h, heterozygotes around 22 h. The mutant gene was later cloned by the group of Joe Takahashi at Northwestern University and found to be equivalent to that coding for a casein kinase enzyme CK1 [START_REF] Lowrey | Positional syntenic cloning and functional characterization of the mammalian circadian mutation tau[END_REF]. Takahashi's group had already discovered the first mammalian clock gene in 1994. This was based upon a mutation, again produced by a mutagenesis screen, causing a long period, and baptized Clock [START_REF] Vitaterna | Mutagenesis and mapping of a mouse gene clock, essential for circadian behaviour[END_REF]. Since then, scientists kept unveiling more details about the transcriptional-translational feedback loops involved in the circadian molecular process. Mathematicians and theoretical biologists pursued the trend and elaborated mathematical models to describe the molecular machinery behind the circadian clock. This ranges from simple negative feedback loop models to complex molecular models [START_REF] Ruoff | Circadian rhythms and protein turnover: the effect of temperature on the period lengths of clock mutants simulated by the goodwin oscillator[END_REF][START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF][START_REF] Gérard | Entrainment of the mammalian cell cycle by the circadian clock: Modeling two coupled cellular rhythms[END_REF][START_REF] Forger | A simpler model of the human circadian pacemaker[END_REF][START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF].

Molecular mechanism of the mammalian circadian clock

This subsection about the molecular mechanism of the circadian clock is based on a book chapter written by E.D. Buhr and J.S. Takahashi [START_REF] Buhr | Molecular components of the mammalian circadian clock[END_REF].

The circadian clock is present in almost all mammalian tissues and there is a central circadian clock in the suprachiasmatic nuclei (SCN) that serves as a "master" clock for the entire body [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF][START_REF] Tosini | Circadian rhythms in cultured mammalian retina[END_REF][START_REF] Brown | Peripheral circadian oscillators in mammals[END_REF][START_REF] Stephan | Circadian rhythms in drinking behavior and locomotor activity of rats are eliminated by hypothalamic lesions[END_REF][START_REF] Moore | Loss of a circadian adrenal corticosterone rhythm following suprachiasmatic lesions in the rat[END_REF][START_REF] Slat | The clock in the brain: Neurons, glia, and networks in daily rhythms[END_REF]. The master clock, through neural, humoral and systemic signals, synchronizes individual cells and tissues to a uninform internal time. This helps maintaining a regular ∼ 24 h rhythm for different processes in mammalian organisms, like body temperature, blood pressure, circulating hormones and metabolism [START_REF] Aschoff | Circadian control of body temperature[END_REF][START_REF] Green | The meter of metabolism[END_REF]. We have seen in the previous section that circadian rhythms are endogenous and encoded in the genetic material. The mechanism behind these rhythms is a coordinated system of transcriptionaltranslational feedback loops that occur within each cell. This system incorporate three main loops. The most essential one is a negative feedback loop in which the heterodimer BMAL1/CLOCK binds to E-boxes and E'-boxes on the promoter region of Per and Cry genes in the nucleus and initiate their transcription. The resulting Per and Cry mRNA migrate to the cytoplasm and translates into proteins. PER and CRY proteins bind together, form a complex and translocate to the nucleus to inhibit BMAL1/CLOCK activity, destroying in this way their transcription factor [START_REF] Griffin | Light-independent role of CRY1 and CRY2 in the mammalian circadian clock[END_REF][START_REF] Sato | Feedback repression is required for mammalian circadian clock function[END_REF]. In addition to that, studies showed that there are two kinases, CK1 and CK1δ that play a role in PER and CRY degradation [START_REF] Lowrey | Positional syntenic cloning and functional characterization of the mammalian circadian mutation tau[END_REF][START_REF] Camacho | Human casein kinase Iδ phosphorylation of human circadian clock proteins period 1 and 2[END_REF][START_REF] Vanselow | Differential effects of PER2 phosphorylation: Molecular basis for the human familial advanced sleep phase syndrome (FASPS)[END_REF]. This kinase-mediated degradation turned out to be of primary importance for terminating the repression phase and restarting the transcription process. For the second feedback loop, BMAL1/CLOCK initiate the transcription of Rev-erbα/β and RORα/β genes [START_REF] Preitner | The orphan nuclear receptor REV-ERBα controls circadian transcription within the positive limb of the mammalian circadian oscillator[END_REF][START_REF] Sato | A functional genomics strategy reveals Rora as a component of the mammalian circadian clock[END_REF][START_REF] Guillaumond | Differential control of Bmal1 circadian transcription by REV-ERB and ROR nuclear receptors[END_REF]. The resulting proteins play different roles, REV-ERB proteins inhibit Bmal1 transcription while ROR proteins initiate it. It was thought that this loop only played a minor role in the sustained circadian oscillations. However, a recent study showed that Rev-erbs are necessary for normal period regulation of circadian rhythmicity [START_REF] Cho | Regulation of circadian behaviour and metabolism by REV-ERB-α and REV-ERB-β[END_REF]. In the third loop, BMAl1/CLOCK binds to D-box elements on the promoter region of a set of PAR bZIP genes, initiating their transcription. This includes genes from the HLF, DBP, TEF and Nfil3 families [START_REF] Falvey | The rat hepatic leukemia factor (HLF) gene encodes two transcriptional activators with distinct circadian rhythms, tissue distributions and target preferences[END_REF][START_REF] Lopez-Molina | The DBP gene is expressed according to a circadian rhythm in the suprachiasmatic nucleus and influences circadian behavior[END_REF][START_REF] Mitsui | Antagonistic role of E4BP4 and PAR proteins in the circadian oscillatory mechanism[END_REF]. Until now, this loop is considered as unnecessary for circadian clock function, but some studies suggest that it helps making the oscillation more robust and precise [START_REF] Liu | Redundant function of REV-ERBα and β and non-essential role for Bmal1 cycling in transcriptional regulation of intracellular circadian rhythms[END_REF][START_REF] Preitner | The orphan nuclear receptor REV-ERBα controls circadian transcription within the positive limb of the mammalian circadian oscillator[END_REF]. Although the transcription-translation feedback loop of Per and Cry genes plays the central role in generating sustained oscillations, it has been speculated that the other feedback loops may provide the necessary delay to create the ∼ 24 h period.

An important distinction between the central clock and peripheral clocks is the difference of the effects of genes mutation on the clock behavior. SCN with Clock knockout, for example, was shown to remain rhythmic [START_REF] Debruyne | A clock shock: Mouse CLOCK is not required for circadian oscillator function[END_REF]. This is because the gene Npas2, plays the role of Clock in its absence as a transcriptional partner for Bmal1 in initiating the transcription activity. However, Clock knockout abolishes rhythmic activities in peripheral clocks [START_REF] Debruyne | Peripheral circadian oscillators require CLOCK[END_REF]. Furthermore, The SCN was shown to stay rhythmic in the absence of any of Per or Cry genes. This turns out to be quite different for peripheral clocks; Cry2 mutation implies a clock with longer period, and Cry1, Per1, Per2 lead to an arrhythmic clock [START_REF] Liu | Intercellular coupling confers robustness against mutations in the SCN circadian clock network[END_REF].

Circadian clock models

The subsection about phase response curve models is adapted from The Colin S. Pittendrigh Lecture written by S. Daan [START_REF] Colin Pittendrigh | and the natural entrainment of circadian systems[END_REF] and the work of J.J. Tyson and L. Glass written for the memorial of Art Winfree [152]. The rest of the subsection is based on our published paper El Cheikh el al. [START_REF] Cheikh | Modeling biological rhythms in cell populations[END_REF].

Phase response curve models

Early works (1950s, 1960s and early 1980s) focused on studying the entrainment properties of the circadian clock with the help of phase response curve models [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF][START_REF] Winfree | Biological rhythms and the behavior of populations of coupled oscillators[END_REF][START_REF] Wever | Zum Mechanismus der biologischen 24-Stunden-Periodik[END_REF][START_REF] Pittendrigh | Circadian systems: Entrainment[END_REF].

A phase response curve or a phase resetting curve is established by delivering a precisely timed perturbation to an oscillation and measuring its effects on the cycle period and phase. The focus, by that time, was on the impact of light delivery on the endogenous circadian oscillations. There were two main theoretical approaches for the nature of delivered light that led to two types of entrainment: parametric and non-parametric. The principle of parametric entrainment, which was initiated by Aschoff, is based on delivering a continuous or tonic light and studying the changes in period as a result of different light intensities [START_REF] Aschoff | Angeborene 24-Stunden-Periodik beim Käcken[END_REF].

The principle of non-parametric entrainment, which was initiated by Colin S. Pittendrigh, is based on delivering a discrete light pulse and studying its effects on the oscillations large phase shifts [START_REF] Pittendrigh | Circadian systems: Entrainment[END_REF]. Pittendrigh thought that the free running period τ of the circadian system was corrected each day for the difference between τ and 24 h when the light fell at a particular phase of the cycle in which a phase shift ( φ = τ -T ) is generated. Aschoff suggested that light lengthens or shortens the period of an endogenous oscillation, while affecting at the same time the average level around which this oscillation moves. He suggested that light may affect the period of the circadian oscillation and at the same time modify its shape and the level around which the oscillation moves. These two different points of view for entrainment emanate from the different fields Pittendrigh and Aschoff were working on.

Pittendrigh elaborated his studies on the single instantaneous event of Drosophila eclosion, whereas Aschoff was interested in modulation of circadian activity of birds, mammals and humans. Following Aschoff work, his collaborator Rütger Wever constructed a mathematical model to study light entrainment of circadian oscillations. His model was inspired from a model originally used for electrical circuits; the van der Pol oscillator [START_REF] Van Der Pol | A theory of the amplitude of free and forced triode vibrations[END_REF]. The van der Pol oscillator is a particularly simple differential equation that can produce a stable limit cycle:

d 2 y dt 2 + (y 2 -1) dy dt + y = 0. (1.2.1)
This equation models an oscillator y(t) with a nonlinear damping coefficient (y 2 -1). The damping term dictates the dynamics of the oscillator. The equation reduces to the harmonic oscillator when = 0. For positive values of , the system is non conservative, and a limit cycle exists. The existence of the limit cycle is verified by noticing that when |y| > 1, the damping is positive (the oscillator dissipates energy) and the amplitude of y decreases. When |y| < 1, the damping is negative (the system receives energy), and the amplitude of y increases. The limit cycle is the trajectory for which the average energy balance is null. Wever model was the following

d 2 S dt 2 + 0.5(S 2 + S -2 ) dS dt + S + 0.6S 2 = L + dL dt + d 2 L dt 2 (1.2.2)
Using this model, Wever studied the effects of light intensity L on the period of oscillations and on their waveform. He concluded that increasing light intensity increases the oscillations frequency as well as their activity time which he denoted by α. He also depicted different waveforms for different photoperiods.

Another example of phase curve models comes from the eminent work of Art Winfree who worked in Pittendrigh lab. Unlike other theoretical biologists who studied the case of a single linear oscillator being forced by an external force, Winfree studied how multiple oscillators can be brought together to oscillate in a synchronized way. He concluded that a weak coupling keeps every oscillator with its autonomous period and the population of oscillators oscillate in an asynchronous fashion. However, he found a threshold coupling strength that made the population of oscillators transit from an asynchronous state to a synchronous one where all oscillators oscillate with the same period [START_REF] Winfree | Biological rhythms and the behavior of populations of coupled oscillators[END_REF]. Like his predecessors, Winfree also studied the impact of light entrainment on circadian oscillations. He plotted the "new phase" (after a light pulse) as a function of the "old phase". He showed that a weak pulse causes only a small change in phase, while a strong pulse drives the oscillator to almost a constant phase regardless of the autonomous phase. Due to the periodicity of the phase variable, he had the brilliant idea of plotting the phase curves on a torus and identified two types of phase resetting curves. Type 1: for a weak pulse, the line relating new phase to old phase passes once through the hole of the torus. Type 0: for a strong pulse, the line does not pass through the hole. Winfree believed that the function relating new phase to old phase and light intensity cannot be continuous; for if it was continuous, the topological invariant of the curve could never jump discontinuously from 1 to 0. This led him to search for a critical value (old phase critical, light intensity) for which the function relating the new phase to these critical values is undefined. He thought that by finding these critical values, he can desynchronize any population of oscillators and he constructed a "fly machine" to test his predictions on the circadian rhythms of fruit flies. By delivering short light pulses at precise duration and intensity, he found the critical values (old phase critical, light intensity) that desynchronize the flies eclosion rhythm [START_REF] Winfree | Integrated view of resetting a circadian clock[END_REF].

Goodwin-type models

Van der Pol type models have been influential in the circadian modeling literature [START_REF] Achermann | Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves[END_REF][START_REF] Czeisler | Bright light induction of strong (type 0) resetting of the human circadian pacemaker[END_REF][START_REF] Forger | A simpler model of the human circadian pacemaker[END_REF][START_REF] Forger | Reconciling mathematical models of biological clocks by averaging on approximate manifolds[END_REF][START_REF] Kunz | Simulation of circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators[END_REF]. However, analysis of single-cell imaging studies from the past decade suggested that cell oscillators could be sloppy, or even damped [START_REF] Yamaguchi | Synchronization of cellular clocks in the suprachiasmatic nucleus[END_REF][START_REF] Westermark | Quantification of circadian rhythms in single cells[END_REF][START_REF] Welsh | Suprachiasmatic nucleus: Cell autonomy and network properties[END_REF][START_REF] Webb | Intrinsic, nondeterministic circadian rhythm generation in identified mammalian neurons[END_REF][START_REF] Maywood | Synchronization and maintenance of timekeeping in suprachiasmatic circadian clock cells by neuropeptidergic signaling[END_REF][START_REF] Bernard | Synchronization-induced rhythmicity of circadian oscillators in the suprachiasmatic nucleus[END_REF]. This is incompatible with the van der Pol model, which always produces limit cycle oscillations (for positive value of ). Theoretical biologists and mathematicians started developing biochemicallybased models. One of the first and probably most popular biochemically-based models was the Goodwin model [START_REF] Goodwin | Temporal Organization in Cells. A Dynamic Theory of Cellular Control[END_REF][START_REF] Goodwin | Oscillatory behavior in enzymatic control processes[END_REF]. The Goodwin model refers to a class of generic molecular oscillators based on a negative feedback loop (the final product of a three-step chain of reactions inhibits the production of the first component) [START_REF] Gonze | Modeling circadian clocks: From equations to oscillations[END_REF]. The original Goodwin equations are

dx dt = k 0 k n 1 k n 1 + z n -k 2 x, (1.2.3) dy dt = k 3 x -k 4 y, (1.2.4) dz dt = k 5 y -k 6 z. (1.2.5)
All equations are linear except for the first one. The nonlinear term in the first equation is a negative feedback term, called a Hill function, and z acts negatively on the production of x. The repressor z can be viewed as a delayed version of the variable x. Standard linear stability analysis shows that three variables in the Goodwin model are necessary for a limit cycle to exist. With a relatively large Hill coefficient (n), the system can oscillate. Thus, in addition to the negative feedback loop, a delay is a necessary ingredient to obtain sustained oscillations. Large Hill coefficients are usually not biologically realistic, but if more intermediate variables or more nonlinear terms are added, the Goodwin model can oscillate for smaller Hill coefficients [START_REF] Gonze | The Goodwin model: Behind the Hill function[END_REF]. Either modifications to the model make it more complex to analyze, and details about intermediate steps are often unknown. Modeling circadian rhythms should take into account all the transcriptional and translational activities that are behind the oscillatory phenomena of the circadian clock, this will lead to a complex set of equations with a lot of parameters.

Instead of detailing all intermediate processes, it is tempting to introduce a "time delay" on x that takes into account the time required to produce the repressor z. This time delay can be introduced in a clean way into the Goodwin model. Let x be the amount of an activator (for example the concentration of mRNA or a protein), which produces through a linear chain process a quantity z, which in turn regulates x. We suppose that the regulator z is the product of a linear chain of differential equations of length p, with kinetic parameter a: ). The gamma density can be viewed as the distribution of time required for the signal activated by x to affect the production of x. This distribution is characterized by a mean delay τ = p/a and a variance p/a 2 . The number of steps p and the kinetic rates thus determine the position and the shape of the delayed distribution. When the number of steps in the linear chain p and the kinetic rates a go to infinity while the mean is constant, the gamma density converges to a Dirac mass at τ , and z(t) = x(tτ ). Hence, feedback loops with large number of intermediate steps can be described with a discrete delay differential equation of the form

dy 1 (t) dt = a x(t) -y 1 (t) , ( 1
dx(t) dt = k 0 k n 1 k n 1 + [x(t -τ )] n -αx(t). (1.2.11)
This formulation of a negative feedback loop is a convenient way to capture the key role of the linear activation chain: producing a delay. This delay is necessary for a limit cycle to exist. Setting τ = 0 in equation ( 1

Molecular models

During the 90s, biologists started to unveil the molecular mechanism of the circadian clock. This helped computational biologists to elaborate more realistic models, based on the molecular processes that generate the sustained circadian oscillation. One of the first models in this category is the model constructed by Albert Goldbeter in 1995 [START_REF] Goldbeter | A model for circadian oscillations in the drosophila period protein (PER)[END_REF]. By that time, it was known that an important part of the mechanism of circadian rhythms is the circadian variation of Per mRNA. The corresponding PER proteins follow a delayed version of Per mRNA rhythms; suggesting that circadian oscillations involve a negative feedback exerted by PER on the transcription of Per gene [START_REF] Hardin | Feedback of the drosophila period gene product on circadian cycling of its messenger RNA levels[END_REF][START_REF] Hardin | Circadian oscillations in period gene mRNA levels are transcriptionally regulated[END_REF]. It was also shown that this delay is caused partly by the post-translational regulation of PER proteins through multiple phosphorylation cascades. Ruoff et al. proposed that this loop could be described by a generic Goodwin type system. They associated the three equations of the Goodwin model to the aforementioned process where a gene mRNA is translated into a protein that activates a transcription factor, which itself inhibits its own gene [START_REF] Ruoff | Circadian rhythms and protein turnover: the effect of temperature on the period lengths of clock mutants simulated by the goodwin oscillator[END_REF].

Goldbeter's model contains more details and is based on the phosphorylation cascades of PER and on the negative feedback exerted by PER on its own gene. We give here some details about his model, because it serves as a prototype for the construction of more complex models including more variables and molecular processes. The model equations are the following:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dM dt = v s K n 1 K n 1 + P n N -v m M K m + M , dP 0 dt = k s M -V 1 P 0 K 1 + P 0 + V 2 P 1 K 2 + P 1 , dP 1 dt = V 1 P 0 K 1 + P 0 -V 2 P 1 K 2 + P 1 -V 3 P 1 K 3 + P 1 + V 4 P 2 K 4 + P 2 , dP 2 dt = V 3 P 1 K 3 + P 1 -V 4 P 2 K 4 + P 2 -k 1 P 2 + k 2 P N -v d P 2 K d + P 2 , dP N dt = k 1 P 2 -k 2 P N .
(1.2.12)

In this model, the variable M is Per mRNA and is considered to be synthesized in the nucleus and transfer to the cytosol where it is degraded. P 0 , P 1 , P 2 are respectively the unphosphorylated, mono-phosphorylated and bi-phosphorylated forms of PER protein. P N is the nuclear form of PER protein. Several assumptions regarding the molecular processes were made: the rate of translation of PER is proportional to M ; P 2 is marked both for degradation and reversible transport into the nucleus; degradation of PER could also be directed at the nuclear form P N as well as to the unphosphorylated or monophosphorylated forms of the protein that could both be transported into the nucleus. This cascade of phosphorylation generate the delay necessary to create sustained oscillations. Also, the first equation is an equation of Hill type describing the negative feedback loop exerted by P N on Per gene and is crucial for the mechanism of oscillations.

This model produces sustained periodic oscillations for a wide range of parameter values. The phase shift between the peaks of total PER proteins and Per mRNA is around 4.5 h. This is in agreement with the experimental results obtained by Zeng et al 1994 [START_REF] Zeng | Constitutive overexpression of the Drosophila period protein inhibits period mRNA cycling[END_REF]. Sustained oscillations in PER and Per mRNA correspond to the evolution toward a stable limit cycle (Figure 1

.3B).

A comprehensive model for the mammalian circadian clock was later developed by Leloup and Goldbeter [START_REF] Leloup | Toward a detailed computational model for the mammalian circadian clock[END_REF]. They have incorporated the regulatory effects exerted on genes expression by the PER, CRY, BMAL1, CLOCK and Rev-ERBα proteins. In a simplified way, the model stated that the oscillations in the core of the circadian system can be generated through two negative feedback loops. The first one exerted on the expression of Per and Cry genes through the Binding of PER-CRY to the CLOCK-BMAL1 activated complex. The second one exerted by CLOCK-BMAL1 through REV-ERBα on the expression of Bmal1 gene. A more detailed model was given by Forger and Peskin who made the distinction between two categories of PER proteins: PER1 and PER2 and the CRY proteins: CRY1 and CRY2. They have used a slower rate of phosphorylation for PER1 because it requires more phosphorylation to bind with CRY1 or CRY2. They also used a higher coefficient of degradation for CRY2 because it is ubiquitinated more quickly than CRY1 [START_REF] Forger | A detailed predictive model of the mammalian circadian clock[END_REF]. Both models have suggested that light can enhance the transcription activity by inducing the production of PER mRNA.

The model that captures most the reality of circadian clock molecular processes is perhaps the one developed by Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF]. Compared to the detailed models of Leloup et Finally, this model is one the best models to fit mutation phenotypes data. It has been observed that dispersed SCN neurons or peripheral tissues/cells lack functional intercellular coupling and thus display independently phased, cell-autonomous rhythms. The intact SCN forms a network of neurons that are synchronized through cell-cell communication mediating intercellular coupling that makes the SCN rhythms robust; but mask cell autonomous circadian phenotypes. Therefore, the effect of knockout mutations on phenotype can be different depending on whether ones assesses at the cellular level or at the SCN tissue or organism level. For example, SCN explants of both Per1 and Cry1 knockouts retain rhythmicity, whereas dispersed SCN neurons of both knockout types are largely arrhythmic [START_REF] Liu | Intercellular coupling confers robustness against mutations in the SCN circadian clock network[END_REF].

In the end of this section, we review in detail a model for the mammalian circadian clock developed by Becker-Weimann et al. [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF] that will play an important role in the coming chapters. It is a good example of a model which takes into consideration real molecular information but has the advantage of being simple, including few variables and equations, and still capturing essential features of the circadian clock functions. The model explores the interdependence of the positive and negative feedback loops created by the transcription factor BMAL1/CLOCK. To make the model less complex, the authors made several simplifying assumptions regarding the molecular processes and lumped multiple circadian actors together. They only considered the activation by BMAL1, relying on the fact that CLOCK is expressed at a constant level, using it as a fixed parameter. CRY proteins are represented by a combined variable, and only Per2 is considered. Per2, Cry mRNAs and proteins are represented by the same variables for several reasons: their expression is coregulated by BMAL1/CLOCK, they form a complex that is necessary for nuclear accumulation, they are both targets of casein kinase I /δ and they both act negatively on BMAL1/CLOCK transcription activity. CKI is considered implicitly by assuming fast phosphorylation of PER2 and CRY. Finally, the repression of Bmal1 by REV-ERBα is taken into account implicitly, by assuming a positive action of PER2/CRY complex on Bmal1 transcription.

The model is given by a set of 7 equations, the first one reads

dy 1 dt = f 1 (y 3 , y 7 ) -k 1d y 1 , (1.2.13)
where:

f 1 (y 3 , y 7 ) = ν 1b (y 7 + c) k 1b (1 + ( y 3 k 1i ) p ) + y 7 + c .
The variable y 1 represents the concentration of P er2 or Cry mRNA which are considered to be identical, for the aforementioned reasons. As one can see from the expression of f 1 (y 3 , y 7 ), the activated form of BMAL1 y 7 , activates the transcription of Per2/Cry mRNA. Whereas increasing the nuclear concentration of PER2/CRY protein (y 3 ), decreases the rate of Per2/Cry mRNA. The coefficient c plays the role of a switch-like behavior of this transcriptional regulation, ν 1b is the maximal rate of transcription, k 1b is the Michaelis constant and k 1d is the degradation rate. The second equation reads

dy 2 dt = k 2b y q 1 -k 2d y 2 -k 2t y 2 + k 3t y 3 , (1.2.14)
The variable y 2 represents the concentration of cytoplasmic PER2/CRY complex. The coefficient k 2b is the rate of formation and k 2d is the rate of degradation of the complex. The coefficients k 2t and k 3t represent respectively the nuclear import and export of PER2/CRY, this justifies the negative sign in front of k 2t and the positive sign in front k 3t . The third equation reads

dy 3 dt = k 2t y 2 -k 3t y 3 -k 3d y 3 , (1.2.15)
Here, the variable y 3 represents the nuclear concentration of PER2/CRY complex. This justifies the opposite signs in front of k 2t and k 3t compared to the previous equation. The coefficient k 3d represents the degradation rate of the complex. The fourth equation reads

dy 4 dt = f 2 (y 3 ) -k 4d y 4 , (1.2.16)
The variable y 4 represents the concentration of Bmal1 mRNA, its rate of transcription is given by:

f 2 (y 3 ) = ν 4b y r 3 k r 4b + y r 3 .
One can see that the transcription rate of Bmal1 increases with rising PER2/CRY (y 3 ) concentration (see Figure 1.3A, dashed line for PER2/CRY and dotted line for Bmal1 ). The coefficient k 4d is a degradation rate. It is noteworthy here to recall that this positive action of PER2/CRY describes the repression process of Bmal1 transcription by REV-ERBα. Hence, this latter protein is included implicitly in the model. The fifth equation reads

dy 5 dt = k 5b y 4 -k 5d y 5 -k 5t y 5 + k 6t y 6 , (1.2.17)
The variable y 5 represents the concentration of cytoplasmic BMAL1 protein. The coefficient k 5b is a translation rate, k 5d is the degradation rate and the coefficients k 5t and k 6t represent respectively the nuclear import and export of BMAL1. The sixth equation reads

dy 6 dt = k 5t y 5 -k 6t y 6 -k 6d y 6 + k 7a y 7 -k 6a y 6 , (1.2.18)
The variable y 6 represents the concentration of nuclear BMAL1 protein. The coefficients k 6d , k 6a and k 7a represent respectively, the degradation rate, the activation and deactivation of BMAL1. The last equation reads

dy 7 dt = k 6a y 6 -k 7a y 7 -k 7d y 7 , (1.2.19)
The variable y 7 represents the activated form of BMAL1 (usually noted BMAL1 * ), which can be understood as its phosphorylated form or its combination with CLOCK. This model shows sustained oscillations for the circadian components with a period approximately equal to 24 hours (Figure 1.3A).

Cell cycle

This section is based on our published paper R. El Cheikh et al. [START_REF] Cheikh | Modeling biological rhythms in cell populations[END_REF]. The cell division cycle is one of life's defining attributes by which an organism is able to reproduce and perpetuate its own species. It can be described as the process where a new born cell doubles its size, replicates its genetic material as well as its cellular components and gives birth to two new progeny cells who inherit all the machinery and information needed to repeat the process.

The cell cycle consists of four distinct phases: G1 (Gap 1), S (DNA synthesis), G2 (Gap 2), and M (mitosis). The G1 phase, is a growth phase, where the cell grows in size and accumulates all the nutrients needed to start DNA synthesis. The S phase is the phase where DNA replication starts, every DNA molecule is replicated into two sister chromatids. When this process is terminated, the cell enters the G2 phase, which acts like a checkpoint. During this gap, the cell ensures that DNA replication has been well accomplished and prepares for entry into mitosis. The mitotic phase consists of four subphases: prophase, where the nuclear envelope is broken and a mitotic spindle is formed; metaphase, where all chromosomes are aligned at the middle of the spindle, with sister chromatids attached to opposite poles of the spindle; anaphase, where the sister chromatids are separated and migrate to the opposite sides of the cell; telophase where two nuclear envelopes are formed and sequester the two complete sets of unreplicated chromosomes. Finally the cell exits mitosis by dividing into two new daughter cells in the G1 phase.

The cell division cycle is a very stringent process of irreversible successions, triggered by transient signals, between its four phases. Progression through the cycle stops if these phases do not take place in the right order. In particular, DNA replication and chromosomes segregation should alternate in proliferating cells. There exist check-point controls that verify whether the steps of the cycle are taking place in the right order. If a problem arises, for example when the replicated chromosomes have not properly aligned on the mitotic spindle, the cell cycle will never exit mitosis. If the chromosomes are well aligned, the check point condition is satisfied, and the transition to the next phase is triggered by transient signals.

These signals disappear once the cell transits to the next phase, making these transitions irreversible.

There are two main types of cell cycle models, molecular models and population models. The first one attempts to model the molecular events of the cell cycle while the second one attempts to describe the dynamics of a population of cells, with an emphasis on cell division (birth) and death events. Before introducing mathematical models, we describe the molecular mechanism of the cell cycle in the coming paragraph.

Molecular mechanism of the cell cycle

In 2001, Paul Nurse, Timothy Hunt and Leland Hartwell won the Nobel Prize for Medicine for their seminal discoveries of key molecular regulators of the cell cycle. They associated the alternation of cell cycle phases to the fluctuation of certain enzyme activity.

Nurse identified specific enzymes called the cyclin-dependent kinases CDKs. These enzymes drive the cell cycle by catalyzing phosphorylation of proteins crucial for the cell cycle progression. Hunt discovered a group of proteins that bind to CDKs to form complexes which are thoroughly required for their activation. At the time, Hunt called newly discovered proteins cyclins, because he liked cycling. The name cyclin turned out to be quite appropriate, since cyclin concentration varies periodically during cell cycle. This implies also the periodic variation of CDKs activity.

CDKs activity is governed, in general, by three distinct mechanisms [START_REF] Morgan | Principles of Cdk regulation[END_REF]. First, regulation is provided by cyclin availability; kinase subunits are present in excess during the cell cycle, but they have no activity until they bind to a cyclin partner. The availability of cyclin subunits is strictly controlled by transcription factors that regulate the expression of cyclin genes, and by ubiquitin-dependent proteolysis systems (e.g., the anaphase-promoting complex APC), which can rapidly degrade cyclin proteins in response to specific signals [START_REF] Peters | The anaphase promoting complex/cyclosome: A machine designed to destroy[END_REF]. Second, CDK activity is regulated by phosphorylation of kinase subunits. Active cyclin/CDK dimers can be inactivated by phosphorylation on a specific tyrosine residue close to the amino terminus of the kinase polypeptide chain. This tyrosine residue is phosphorylated by kinases of the Wee1 family and dephosphorylated by phosphatases of the Cdc25 family [START_REF] Coleman | Cdc2 regulatory factors[END_REF]. Finally, active cyclin/CDK dimers can also be inactivated by binding to inhibitors, called CDK inhibitors (CKIs) [START_REF] Donovan | P40SDB25, a putative CDK inhibitor, has a role in the M/G1 transition in Saccharomyces cerevisiae[END_REF][START_REF] Moreno | Regulation of progression through the G1 phase of the cell cycle by the rum1+ gene[END_REF][START_REF] Schwob | The B-type cyclin kinase inhibitor p40SIC controls the G1 to S transition in S. cerevisiae[END_REF]. The levels of CKIs depend on their production rate, which is governed by regulated transcription factors and their destruction rate (phosphorylated CKIs are rapidly ubiquitinated and degraded) [START_REF] Benito | Regulation of the G1 phase of the cell cycle by periodic stabilization and degradation of the p25rum1 CDK inhibitor[END_REF][START_REF] Verma | Phosphorylation of Sic1p by G1 Cdk required for its degradation and entry into S phase[END_REF].

In molecular models, proteins and mRNA concentrations are often modeled with ordinary differential equations. Models are built with weakly connected modules, each module standing for one of the major checkpoints of the cell cycle. Conceptually, the cell cycle stops at each checkpoint and progression is halted until all conditions are met to raise the checkpoint. The cell cycle can be divided into three main checkpoints: START, G2/M and EXIT. The START checkpoint ensures that the cell has grown enough and have all the sufficient nutrients to leave G1 phase and dedicate to a new round of DNA synthesis. The G2/M checkpoint ensures that there is no damaged DNA after replication before cell enters mitosis. The Exit checkpoint ensures that the chromosomes are properly aligned and that the cell is ready to divide into two healthy daughter cells.

The progression through these checkpoints is triggered by CDK activity fluctuation. In G1 phase, CDK activity is low. To pass the START checkpoint, CDK activity rises up to a certain threshold that allows the cell to start DNA replication process and enter mitosis. To pass the EXIT checkpoint, CDK activity must fall to let the cell divides and the daughter cells to enter G1.

So what causes this fluctuation of CDK activity? There exist protein enemies that block CDK activity. This makes each checkpoint like a battle between CDK and its enemies. During the START transition, CDK activity is downregulated by the enemies and; hence the cell is blocked in G1 until some helpers called Starter kinases (SK) come for the rescue of CDK by downregulating CDK enemies allowing the cell to pass the START checkpoint. During the EXIT checkpoint, Exit proteins (EP) upregulate CDK's enemies, promoting the CDK abrupt degradation and cell division.

Experimental and theoretical biologists have shown that checkpoints are marked by abrupt transitions, or switches, during which specific cell cycle proteins get quickly activated or deactivated [START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell[END_REF][START_REF] Chauhan | Regulation of mammalian cell cycle progression in the regenrating liver[END_REF][START_REF] Novak | Mathematical model of the cell division cycle of fission yeast[END_REF][START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF][START_REF] Ferrell | Modeling the cell cycle: Why do certain circuits oscillate?[END_REF]. A switch-like transition is a response to a change in concentration of a stimulus that affects the system. When the concentration of the stimulus reaches up a certain threshold, the system switches from one state to another. We distinguish two types of switch-like transition: continuous (sigmoidal switches) and discontinuous (bistable switches). The response in the sigmoidal switch is graded and reversible. By graded we mean that the response increases continuously with signal strength. By reversible we mean that if the signal S increases, reaches its threshold value and the system switches on to the next state, a decreasing value for S can bring back the system to its off state.

The discontinuous responses can further be separated into two kinds: the one-way switch (figure 1.4 (A)) and the toggle switch (figure 1.4 (B)). For the one-way switch, there exists a critical value (S crit ) that the signal should reach to let the response passes to the upper state. Now, if S decreases the response does not fall back to the lower state and stays high, this is why the switch is called in this case irreversible. Notice that for S between 0 and S crit , the system has two stable steady-state responses (lower and upper state) separated by an unstable one. This is why we call the switch bistable. A good example of this irreversible (A) [START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell[END_REF] (see Box 1, Figures 1.e and 1.f in the reference for more details). (A) One way switch: the signal should attain the critical value "Scrit" so that the response passes to the upper state. If the signal decreases, the response stays on the upper state and does not fall back to the lower one. (B) Toggle switch: unlike the one way switch, if the response is on the upper state and the signal decreases, the response falls back to the lower state. switch is apoptosis [START_REF] Legewie | Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability[END_REF]. For the toggle switch, if the response is on the upper state, and the value of S decreases enough, the switch will go back to the lower state. For intermediate values of S (S crit1 < S < S crit2 ), the response can be on the upper or lower state, depending on how S was changed. This sort of two-way discontinuous switch is also called hysteresis. A good example of hysteresis response is the activation of the mitosis promoting factor MPF, or the START and FINISH transitions that will be explained in the coming paragraph.

By coupling sequentially many of those bistable switches, it is possible to devise models that can follow the progression of the cell cycle [START_REF] Tyson | Temporal organization of the cell cycle[END_REF][START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF][START_REF] Chauhan | Regulation of mammalian cell cycle progression in the regenrating liver[END_REF].

Molecular models for the cell cycle

Goldbeter minimal model for mitosis

Albert Goldbeter constructed a minimal model for cell cycle mitotic oscillator in early amphibian embryos [START_REF] Goldbeter | A minimal cascade model for the mitotic oscillator involving cyclin[END_REF]. It is based on the cascade of phosphorylation-dephosphorylation cycles involving cyclin and cdc2 kinase. The system of equations contains only three variables which are cyclin, active cdc2 and active cyclin protease. The basic assumption is that cyclin is synthesized at a constant rate and triggers the activation of cdc2 kinase. The latter activates the cyclin protease which elicits cyclin degradation. This switches off the activation of cdc2 creating a negative feedback loop and resetting the system for a new mitotic cycle.

The system of equations read:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dC dt = v i -v d X C k d + C -K d C dM dt = V 1 (1 -M ) K 1 + (1 -M ) -V 2 M K 2 + M dX dt = V 3 (1 -X) K 3 + (1 -X) -V 4 X K 4 + X with V 1 = C Kc+C V M 1 , V 3 = MV M 3 .
(1.3.20)

In the above equations, C denotes the cyclin concentration, while M and X represent the fraction of active cdc2 kinase and the fraction of active cyclin protease; (1 -M ) thus represents the fraction of inactive (i.e., phosphorylated) cdc2 kinase, while (1-X) represents the fraction of inactive (i.e., dephosphorylated) cyclin protease.

The main result of this model is the demonstration the negative feedback provided by cdc2-induced cyclin degradation suffices to create sustained oscillations with correct period and amplitude. No positive feedback is needed, sustained oscillations arise from the thresholds and time delays built into the cascade of post translational modification controlling the activation of cdc2 kinase and cyclin proteolysis.

M T (total amount of cdc2 kinase) for enzymes E 1 and E 2 , and X T (total amount of cyclin protease) for enzymes E 3 and E 4 ; both M T [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF][START_REF] Bernard | How to build a multiscale model in biology[END_REF][START_REF] Bernard | Tumor growth rate determines the timing of optimal chronomodulated treatment schedules[END_REF] and X T will be considered as constant throughout the cell cycle. The expressions for the effective maximum rates V 1 and V 3 are given by Eq. 2. These expressions reflect the assumption that cyclin activates phosphatase E 1 according to a Michaelis-Menten kinetic; V M 1 denotes the maximum rate of that enzyme reached at saturating cyclin levels. On the other hand, the effective maximum rate of cdc2 kinase is proportional to the fraction of active enzyme; V M 3 denotes the maximum velocity of the kinase reached for M = 1.

Tyson and Novak fission yeast model

Tyson and Novak fission yeast model underlines the main molecular events behind the progress of the cell cycle. The core of this model is based on the activity of the cyclindependent protein kinases complexes Cdc2/Cdc13 (also called MPF or mitosis promoting factor), which are the engine needed to start DNA replication and mitosis. In this model, the cell cycle is punctuated by three transitions: Start, G2/M and Finish. These transitions depend on the concentration of Cdc2/Cdc13 and their enemies. If the activity of Cdc2/Cdc13 is high, the cell progresses through the cell cycle; if it is low, the cell blocks its progression. Each phase transition of the cycle is regulated by specific enemies and helpers that decide whether Cdc2/Cdc13 will win or lose. The Start transition (G1 to S) is governed by the antagonistic interaction between Cdc2/Cdc13 and their enemies Ste9 and Rum1. Ste9 targets Cdc13 to the APC core and promotes their degradation, while Rum1 binds to Cdc2/Cdc13 complexes and inhibits their activity. On the other hand, Cdc2/Cdc13 can also downregulate, by phosphorylation, the activity of Ste9 and Rum1. So what shifts the balance to Cdc2/Cdc13 so that they can win and let the cell passes to the next phase? For the Start transition, there exists starter kinases that help MPF to get the upper hand and phosphorylate Ste9 and Rum1. For the Finish transition, the MPF activity should shut down to let the cell exit mitosis and enter the G1 phase. The helper molecule for this transition is the Slp1/APC complex, which promotes the degradation of Cdc13 and activates Ste9. Hence, the activity of the enemies will win over the activity of MPF, which shuts down and lets the cell exit mitosis. In the G2/M transition, the enemy of MPF is the tyrosine kinase WEE1, which can inactivate Cdc2. To shift the balance toward MPF, a specific phosphatase called Cdc25 removes the inhibitory effect of WEE1. Cdc25 is activated in a positive feedback by MPF (Figure 1.6). The mass M is divided by two in the EXIT phase.

The second equation describes the rate of change in Cdc13/Cdc2 (named Cdc13 T ) complex concentration :

d[Cdc13 T ] dt = k 1 M -k 2 [Cdc13 T ] -k 2 [Ste9][Cdc13 T ] -k 2 [Slp1][Cdc13 T ].
The first term on the right hand side assumes that the rate is proportional to the cell mass, the last three terms are the nonspecific degradation, Ste9 and Slp1 -mediated degradation rates.

The third equation represents the antagonism between WEE1 and MPF through the form of the factor k wee (See auxiliary equations):

d[P reMP F ] dt = k wee ([Cdc13 T ] -[P reMP F ]) -k 25 [P reMP F ] -(k 2 + k 2 [Ste9] +k 2 [Slp1])[P reMP F ],
where PreMPF refers to the activated form of Cdc13 T /Cdc2. The fourth equation describes the rate of change in Slp1 total concentration:

d[Slp1] dt = k 5 + k 5 [MP F ] 4 J 4 5 + [MP F ] 4 -k 6 [Slp1].
The first term on the right hand side is a synthesis term, the second term is a Hill type synthesis term due to MPF, and the last term is a degradation term.

The fifth equation describes the rate of change in the IEP enzyme activity. IEP provides the delay necessary for the chromosomes to align with the metaphase plane before they are separated at anaphase: The sixth equation describes the rate of change of Ste9 total concentration:

d[IEP ] dt = k 9 [MP F ](1 -[IEP ]) J 9 + (1 -[IEP ]) - k 10 [IEP ] J 10 + IEP . (A)
d[Ste9] dt = k 3 + k 3 [Slp1](1 -[Ste9]) J 3 + (1 -[Ste9]) - (k 4[SK] + k 4 [MP F ])[Ste9] J 4 + [Ste9] .
The first term on the right hand side is an activation term and the second one represents deactivation caused by SK and MPF.

The seventh equation describes the variation of the total concentration of Rum1 :

d[Rum1 T ] dt = k 11 -(k 12 + k 12 [SK] + k 12 [MP F ])[Rum1 T ].
The first term on the right hand side is a pure synthesis one, the second term represents constant degradation and degradation due to SK and MPF activities.

The eighth equation represents the variation of the SK kinases concentration:

d[SK] dt = k 13 [T F ] -k 14 [SK]
here TF is some function of the mass M and on MPF.

Renewal equations and structured division models

In this section, we present structured partial differential equations in the context of cell proliferation or population models. This section is based on our published paper R. El Cheikh et al. [START_REF] Cheikh | Modeling biological rhythms in cell populations[END_REF]. We emphasize results related to circadian control, that is, theorems on growth rate with periodic transition coefficients.

Our typical population model is written as it follows

⎧ ⎪ ⎨ ⎪ ⎩ ∂ t n(t, x) + ∂ x n(t, x) = -d(t, x)n(t, x), x > 0, n(t, x = 0) = ∞ 0 B(t, x)n(t, x)dx, ( 1
⎧ ⎪ ⎨ ⎪ ⎩ ∂ t n(t, x) + ∂ x n(t, x) = -[d(t, x) + K(t, x)]n(t, x), x > 0, n(t, x = 0) = 2 ∞ 0 K(t, x)n(t, x)dx.
(

In this model, the loss rate is divided into two parts: a death rate d and a division rate K.

The loss of one cell due to the term K is compensated by the creation of two new cells of age x = 0. The age x is not a physiological maturity but is the chronological age since birth.

Using the method of characteristics, one can show that

n(t + x, x) = n(t, 0) exp - x 0 d(t + s, s)ds .
The latter equation is the key to the link between structured population models and delay differential equations. Using the boundary condition, we can derive a Volterra-like integral equation satisfied by n(t, 0)

n(t, 0) = ∞ 0 B(t, x)n(t -x, 0) exp - x 0 d(t + s, s)ds . (1.4.23)
Well posedness for those models is classical in the case of bounded coefficients and initial conditions in L 1 (R + ) [START_REF] Perthame | Transport equations in biology[END_REF].

More detailed cell cycle models can be devised. In [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF] Lepoutre and coworkers considered a cell cycle model with I phases

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t n i (t, x) + ∂ x n i (t, x) + [d i (t, x) + K i→i+1 (t, x)n i (t, x) = 0, 1 ≤ i ≤ I n i+1 (t, 0) = ∞ 0 K i→i+1 (t, x)n i (t, x)dx, n 1 (t, 0) = 2 ∞ 0 K I→1 (t, x)n I (t, x)dx.
(1.4.24) usually I = 4, and I + 1 = 1. Similar to division models, the loss rate d i + K i in the i-th equation contains two terms: a death rate d i and a transition rate K i from phase i to i + 1, which turns out to be a division rate for i = I. This system was studied with the aim to show that tumor growth is enhanced by circadian clock disruption. It was proved, using a convexity result for the dominant eigenvalue (to be introduced below) of the system, that for a disrupted circadian rhythm (averaged coefficients), the dominant eigenvalue is smaller than the cases of controlled circadian rhythm (periodic coefficients). This result would imply that periodic population grows faster. This lead the authors to conclude that disrupted circadian rhythms do not enhance tumor growth directly but rather damages the healthy tissues that fight against it [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF].

Another important class of population models in cell cycle representation are delay differential models. These models are linked to previous partial differential equations. For instance, in [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF], the following discrete delay equation, which was studied in [START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF][START_REF] Mackey | Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis[END_REF], was re-derived from a division model

dp(t) dt = -[d(t) + K(t)]p(t) + 2σ(t)K(t -τ )p(t -τ ). (1.4.25)
Based on age-structured PDEs, more general delay models can be constructed. For example, systems where the discrete delay is replaced by a distribution of delays can be derived rigorously from division models:

dx(t) dt = -a(t)x(t) + ∞ 0 x(t -u)b(t, u)g(u)du. (1.4.26)
The coefficients a(t), b(t, u) are directly related to the coefficients of the division model with I = 2 in the following way:

a(t) = d 1 (t) + K 1 (t), (1.4.27) g(u) = K 2 (u)e -u 0 K 2 (s)ds , (1.4.28) b(t, u) = 2K 1 (t -u)e -t t-u d 2 (s,s-t+u)ds . (1.4.29)
In the coming paragraph we summarize various properties for the growth rate of such equations. Due to the close link between PDE-based and DDE-based models, results are summarized for PDE models only. We choose to write all the theorems on (1.4.21) but mutatis mutandis the results remain true for division models.

Asymptotic behavior with or without periodic forcing

The theoretical results presented hereafter concerning the analysis of asymptotic behavior have been obtained by T. Lepoutre in his thesis [START_REF] Lepoutre | Analysis and modelling of growth and motion phenomenon from biology[END_REF] and have been published [START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF][START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF][START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF] .

The quantity of interest is the growth exponent. We assume that coefficients satisfy conditions ensuring net growth of the population and that the asymptotic behavior is governed by the principle eigenvalue and its associated eigenvector (examples of such conditions can be found in [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF]). For simplicity, we present the results for the renewal equation (1.4.21). We assume there exists a triple (N, Φ, λ), such that

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t N (t, x) + ∂ x N (t, x) = -(d(t, x) + λ)N (t, x), N (t, x = 0) = ∞ 0 B(t, x)N (t, x)dx, -∂ t Φ(t, x) -∂ x Φ(t, x) + (d(t, x) + λ)Φ(t, x) = B(t, x)Φ(t, 0). λ is a positive number, N, Φ are positive functions ∀t, N (t + T, x) = N (t, x), Φ(t + T, x) = Φ(t, x), 1 T T 0 ∞ 0 N (t, x)dx = ∞ 0 N (t, x)Φ(t, x)dx = 1.
The last condition is just a renormalization to ensure uniqueness. Classically, such a λ is unique and governs the growth (detailed results on general relative entropy theory can be found in [START_REF] Perthame | Transport equations in biology[END_REF]). One has, in particular,

d dt ∞ 0 n(t, x)e -λt φ(t, x)dx = 0, d dt ∞ 0 |n(t, x)e -λt -ρ 0 N (t, x)|φ(t, x)dx ≤ 0
where ρ 0 = ∞ 0 n 0 (x)Φ(0, x)dx > 0 is determined only by n 0 . Conditions for existence are not the aim of this section and for this purpose the reader could refer to [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF][START_REF] Lepoutre | Analysis and modelling of growth and motion phenomenon from biology[END_REF].

To analyze the impact of periodic forcing, it is natural to compare the behavior of the periodic system with a time-constant system. The first result in that direction can be summarized as Theorem 1. [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF] Assume in (1.4.21) that the birth rate B does not depend on time. Let λ per be the growth rate of the system (1.4.21) and λ s be the system where d(t, x) has been replaced by its arithmetical average over time

d s (x) = 1 T T 0 d(t, x)dt. (1.4.30)
then the following inequality holds true

λ per ≥ λ s .
A more general result was then obtained concerning the role of the birth rate B [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF]. Surprisingly, it leads to the apparition of a geometrical average of this rate

B g (x) = exp 1 T T 0 log B(t, x)dt . (1.4.31)
The result can be summarized as it follows Theorem 2. [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF] Define λ per as the growth rate of (1.4.21) and λ g the growth rate of (1.4.21) with d replaced by d s and B replaced by B g , then the following inequality holds true

λ per ≥ λ g .
Introducing different types of means can be puzzling but might be better understood by looking closely at the system structure. In equation (1.4.23), we can see that log B(t, x) and x 0 d(tx + s, s)ds have a similar role. This mathematical result cannot be readily exploited for division or cell cycle models, at least concerning division and transition coefficients, since division rate K also appears as a loss rate and a birth rate in (1.4.22). In this case, the averaged model cannot not be characterized as a division model because K would be 50 We tested this result and compared the growth of a cell population with and without periodic forcing. Simulations were coherent with theoretical results, in the sense that there is no general inequality between the growth with and without periodic coefficients (Figure 1.8).

The first two theorems above can be generalized to a convexity property: Theorem 4. [START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF] The growth rate of equation ( 1

d θ = θd 1 + (1 -θ)d 2 , (1.4.32) B θ = (B 1 ) θ (B 2 ) 1-θ . (1.4.33)
Then, from Theorem 4, the following inequality holds:

λ θ per ≥ θλ 1 per + (1 -θ)λ 2 per .
Theorem 4 is a continuous version of the "Jensen"-based Theorem 2. It can also be seen as an extension to periodic systems of an inequality given by Kingman on spectral radius of nonnegative matrices [START_REF] Kingman | A convexity property of positive matrices[END_REF]. This theorem leads to a theoretical justification of chronotherapy. Indeed, consider a drug, given every day at the same time, which side effects on healthy tissue are only represented by an additional death rate d drug (tt adm , a) (the parameter t adm representing the effect of the drug administration time). Equation (1.4.21) is then replaced by

⎧ ⎨ ⎩ ∂ t n + ∂ x n + [d(t, x) + d drug (t -t adm , x)]n(t, x) = 0, n(t, x = 0) = ∞ 0 B(t, x)n(t, x
)dx, to which growth rate λ(t adm ) is naturally associated. In the case of a continuous treatment, the drug induced death rate would be

d cont (a) = 1 T T 0 d drug (t, a)dt = 1 T T 0 d drug (t -t adm , a)dt adm .
Again a growth rate λ cont is associated. As a consequence of the convexity of the growth rate with respect to death rates, we have

1 T T 0 λ(t adm )dt adm ≥ λ cont .
which implies that an average periodic drug delivery is less toxic than a constant delivery.

Chapter 2

Modeling circadian clock-cell cycle interaction effects on population growth rate

Introduction

The work presented in this chapter is based on our published paper R. El Cheikh et al. [START_REF] Cheikh | Modeling circadian clock-cell cycle interaction effects on cell population growth rates[END_REF].

In the previous chapter, we presented the cell cycle and the circadian clock as two distinct biological oscillators, we discussed their molecular mechanisms and mathematical models employed for their study. However, the circadian clock interacts with the cell cycle through multiple molecular pathways [START_REF] Edmunds | Cellular and molecular bases of biological clocks. Models and mechanisms for circadian timekeeping[END_REF][START_REF] Johnson | Circadian clocks and cell division. What's the pacemaker?[END_REF][START_REF] Yang | Circadian gating of the cell cycle revealed in single cyanobacterial cells[END_REF][START_REF] Pando | Coupling cellular oscillators-circadian and cell division cycles in cyanobacteria[END_REF][START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF][START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomus and self-sustained oscillators pass time to daughter cells[END_REF]. The aim of this chapter is to construct a mathematical model that studies the influence of this coupling on the growth of a cell population. Our work is motivated by several epidemiological studies that have shed light on the fact that individuals with disrupted circadian rhythms have increased risk of developing tumorigenic diseases [START_REF] Kubo | Prospective cohort study of the risk of prostate cancer among rotating-shift workers: Findings from the Japan collaborative cohort study[END_REF][START_REF] Filipski | Circadian disruption accelerates liver carcinogenesis in mice[END_REF][START_REF] Filipski | Host circadian clock as a control point in tumor progression[END_REF][START_REF] Filipski | Effects of chronic jet lag on tumor progression in mice[END_REF][START_REF] Fu | The circadian gene Period2 plays an important role in tumor suppression and DNA damage response in vivo[END_REF][START_REF] Gery | The circadian gene Per1 plays an important role in cell growth and DNA damage control in human cancer cells[END_REF][START_REF] Grèchez-Cassiau | The circadian clock component BMAL1 is a critical regulator of p21 WAF1/CIP1 expression and hepatocyte proliferation[END_REF]. Studies made on yeast revealed restriction of cell division to the reductive phases of the yeast metabolic cycle. This regulation insures that the cell cycle evades the potentially mutagenic redox environment of the oxidative respiratory phase, helping to minimize the occurrence of futile reactions [START_REF] Murray | Regulation of yeast oscillatory dynamics[END_REF][START_REF] Tu | Logic of the yeast metabolic cycle: Temporal compartmentalization of cellular processes[END_REF]. This type of control was shown to be involved in circadian regulation and may be a general strategy for the robust maintenance of cellular processes. We have seen, as described in Section 1.3.1, progression through each phase of the cell cycle depends on the activity of cyclins and cyclin-dependent protein kinase complexes (Cdks) and a mitosis promoting factor (MPF). When MPF activity is high, the cell progresses through the cycle. When it is low, progression stops [START_REF] Morgan | Principles of Cdk regulation[END_REF]. Each phase of the cell cycle is controlled by a different cyclin/Cdk complex: G1 is controlled by cyclin D/Cdk4-6, G1/S transition by cyclin E/Cdk2, S phase by cyclin A/Cdk2 and G2/M transition by cyclin B/Cdk1 [START_REF] Hunt | Riding tandem: Circadian clocks and the cell cycle[END_REF]. The link between the circadian clock and the cell cycle is due to the regulation of different Cdks by several molecular clock components. It has been reported that BMAL1/CLOCK activates the transcription of the kinase WEE1 to regulate the G2/M transition [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF]. The circadian clock, via the protein REV-ERBα, regulates the transcription of p21, which inhibits Cdk2 and blocks the G1/S transition [START_REF] Grèchez-Cassiau | The circadian clock component BMAL1 is a critical regulator of p21 WAF1/CIP1 expression and hepatocyte proliferation[END_REF]. The circadian clock is also involved in direct control of DNA damage and apoptosis pathways by virtue of its regulation of Chk2 and other related factors [START_REF] Gery | The circadian gene Per1 plays an important role in cell growth and DNA damage control in human cancer cells[END_REF][START_REF] Chen | A conserved DNA damage response pathway responsible for coupling the cell division cycle to the circadian and metabolic cycles[END_REF].

Two main approaches have been used to model the coupling between the cell cycle and the circadian clock. The first approach is to model the molecular machinery of the cell. It is usually based on ordinary differential equations, where the variables describe the intracellular molecular concentrations of both oscillators. Chauhan and colleagues constructed such a model to account for the regulation of mammalian cell cycle progression and its gating by the circadian clock in the regenerating liver [START_REF] Chauhan | Regulation of mammalian cell cycle progression in the regenrating liver[END_REF]. Zamborszky and colleagues used a minimal model for circadian rhythms coupled to a cell cycle model that had been originally developed for the yeast cell cycle. Their model revealed quantized cell cycles and they suggested that cell size control is influenced by the clock [START_REF] Zamborszky | Computational analysis of mammalian cell division gated by a circadian clock: Quantized cell cycles and cell size[END_REF]. More recently, Gérard and colleagues used a detailed computational model for the Cdk network driving the mammalian cell cycle to study the effect of multiple molecular links to the circadian clock [START_REF] Gérard | Entrainment of the mammalian cell cycle by the circadian clock: Modeling two coupled cellular rhythms[END_REF]. They characterized the domains of autonomous periods where the cell cycle can be brought to oscillate to 24 or 48 h periods, and determined conditions for switching between these two patterns of entrainment.

The second approach is to model a cell population, leaving aside molecular details. This approach is based on PDEs, especially the category of physiologically-structured models, or on individual-based models and cellular automata. In these models, the cell cycle is divided into multiple, discrete phases and the circadian clock is coupled via time-periodic parameters, such as the transition coefficients or phases duration. Altinok and colleagues used a cellular automaton model to examine the entrainment of the cell cycle by the circadian clock [START_REF] Altinok | An automaton model for the cell cycle[END_REF]. Clairambault and colleagues used an age-structured PDE system to model a population of cells under the control of the circadian clock [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF][START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF]. The circadian clock was taken into account through periodic cell cycle phase transition coefficients into the equations (see section 1.4 for more details).

Compared to population models, molecular models capture more details of the fine regulation of the cell cycle, and in particular, can predict the effect of mutations on the cell cycle regulation. However, molecular models rarely describe explicitly dividing cell populations and it is not clear how growth rates are affected by disruptions at the molecular level.

Here, we present a mathematical model that combines the molecular and the population levels, to study the influence of the circadian clock on the growth of a cell population. We study the influence of circadian clock gene mutations on the net growth rate of a dividing population. We show that disruption of circadian rhythms can lead to abnormal proliferation. Depending on autonomous cell cycle properties and the nature of the disruption, circadian clock gene mutations can lead to faster or slower growth rates. We characterize the effect of circadian clock gene mutations, and show that combined molecular/population model brings to the dynamics of cell proliferation a picture more complete than a molecular or population model alone.

Presentation of the model

Our study focuses on the coupling between the cell cycle and the circadian clock through the protein WEE1. The combined molecular/population model consists of two coupled systems of equations: one system of ordinary differential equations that describes the molecular dynamics of the cell cycle and the circadian clock, and one system of partial differential equations that describes the growth of a cell population. The molecular model itself is a coupled system of two core networks, one for the circadian clock, and one for the cell cycle.

According to Nagoshi et al. [START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomus and self-sustained oscillators pass time to daughter cells[END_REF], cultured fibroblasts harbor self-sustained and cell autonomous circadian clocks similar to those operative in the neurons of the suprachiasmatic nuclei. Circadian gene expression continues during cell division and daughter cells resume the rhythms of mother cells after mitosis [START_REF] Murray | Regulation of yeast oscillatory dynamics[END_REF]. A recent computational study examined the effect of cell division of genetic oscillators. It was shown that oscillations are quite resilient to cell division, and that cell-cell heterogeneity appear to be the main source of variability observed experimentally [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF]. Based on these assumption, we assume that cell cycle divisions do not alter the molecular concentration of the circadian components, neither their rhythms.

For the circadian clock, we used the model proposed by Becker-Weimann and colleagues [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF]. It consists of seven nonlinear ordinary differential equations describing the concentrations of Per /Cry mRNA and PER/CRY protein complexes, and Bmal1 /Clock mRNA and protein complexes (y i , i = 1, ..., 7). We recall that this is relatively a simple model that takes into consideration molecular information and that was used to analyze the roles of feedback loops on the oscillatory dynamics (more details are included in Section 1.2.3). This model was used to explore the role of the negative feedback loop created by the transcription factor complex BMAL1/CLOCK that activates the Period and Chryptochrome genes (Per1, Per2, Cry1 and Cry2 ) (Figure 2.1A). After several hours, PER and CRY proteins form a complex in the cytoplasm, go back to the nucleus and downregulate their own synthesis by inhibiting BMAL1/CLOCK. Once the latter protein complex is inhibited, transcription of PER and CRY stops. Hence, BMAL1/CLOCK is no longer inhibited and the cycle starts its process again. The model also includes a positive feedback loop where Bmal1 transcription is positively regulated by PERs and CRYs because the complex PER/CRY also inhibits the transcription of Rev-erbα, which inhibits the transcription of Bmal1.

For the cell cycle, we used a system of three ordinary differential equations based on MPF activity (z i , i = 8, ..., 10). This model was inspired by the model by Tyson and Novak [START_REF] Novak | Mathematical model of the cell division cycle of fission yeast[END_REF] (more details about the model are presented in Section 1.3.2). The core of Tyson and Novak model is based on the activity of the cyclin-dependent protein kinase complexes CyclinB/Cdk1 (also called MPF for Mitosis Promoting Factor), which are the engine needed to start DNA replication and mitosis. The cell cycle is divided into three phases: G1, S/G2, and M. Transitions from one phase to the other depend on the concentration of MPF and its enemies. When the activity of MPF is high, the cell progresses through the cell cycle; when it is low, the cell blocks its progression. Each phase transition of the cycle is regulated by specific enemies and helpers, which decide whether MPF will win or lose. Transition from G1 to S is governed by the antagonistic interaction between MPF and its enemies APC G 1 and CKI. In the G2/M transition, the enemy of MPF is the tyrosine kinase WEE1, which can inactivate Cdk1. At cell division, or M to G1 transition, MPF activity shuts down to let the cell exit mitosis and enter the G1 phase. The helper molecule for this transition is the APC M complex, which promotes the degradation of CyclinB. In the model, three players are included explicitly: MPF (z 8 ), WEE1 (z 9 ) and the inhibitor of MPF (z 10 ). We supposed that cells enter S/G2 phase when MPF increases above a fixed threshold (θ 1 ), enter mitosis (M phase) when MPF activity rises above that of WEE1, and divide when MPF reaches back a low threshold level (θ 2 ), as it happens during mitosis. Even though the cell cycle model presented here is not quantitative due to the small number of kinetic parameters, it still reproduces a correct qualitative behavior of the cell cycle dynamics. Since we were interested in the effects of coupling the circadian clock to the cell cycle through the protein 56 (A) (B) Figure 2.1: (A) Scheme of the circadian clock network: the activated heterodimer BMAL1/CLOCK (BMAL1 , y 7 ) activates Per2 and Cry genes, which produce Per2/Cry mRNA (y 1 ). PER2 and CRY proteins are synthesized and bind in the cytosol to form a complex (y 2 ) to be transported into the nucleus (y 3 ). This complex inhibits the activity of BMAL1/CLOCK complex, thus destroying its own source of transcription and closing the negative feedback loop. The nuclear complex PER2/CRY (y 3 ) also activates Bmal1 transcription, which produces an increase in Bmal1 mRNA (y 4 ), and cytosolic protein concentration (y 5 ). The BMAL1/CLOCK complex is then transported to the nucleus (y 6 ), where it is activated. The activated BMAL1/CLOCK complex (BMAL1 , y 7 ) restarts the activation process of Per2/Cry. (B) Schematic representation of the coupling between the cell cycle and the circadian clock through the protein WEE1 (z 9 ) [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF]. In the model, the cell cycle is divided into three successive phases G1, S/G2, M. Transitions from one phase to another depend on the activity of MPF (z 8 ). For cells to leave G1 and enter S/G2, MPF activity must exceed a fixed threshold θ 1 = 0.09. For cells to leave S/G2 and enter M phase, MPF activity must exceed that of WEE1 (z 9 ). In the M phase, MPF activates its inhibitor (z 10 ), which represses MPF activity, letting it shut down and forcing the cell to exit mitosis. Division occurs once MPF activity reaches a low threshold level θ 2 = 0.06. Coupling between the cell cycle and the circadian clock is achieved by the transcriptional activation of Wee1, which induces WEE1 activity (z 9 ) by the active BMAL1/CLOCK complex (y 7 ). WEE1, we only considered the antagonistic relation of WEE1 and MPF to avoid simulation artefacts, which may come from other interactions that are not related to our study. Our model reproduces well the evolution of MPF activity, which oscillates in an antagonistic way with the activity of WEE1. Once MPF activity surpasses WEE1 activity, it activates its inhibitor (variable z 10 in our model, which can be associated to APC M in the Tyson and Novak model, Figure 1.6 in Section 1.3.1) to help shutting down its own activity and forces the cell to exit mitosis.

The full, 10-variable molecular system, reads For the cell population system, we used age-structured equations as described in Section 1.4. The molecular model entrains the cell population system through cell cycle phase transition rates, which depend on an average molecular state of the cells. The age-structured system tracks the time elapsed by cells in each cell cycle phase [START_REF] Clairambault | Circadian rhythm and tumor growth[END_REF][START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF][START_REF] Doumic | Analysis of a population model structured by the cells molecular content[END_REF]. We divided the cell cycle model into three phases corresponding to the three phases of the molecular model: G1, S/G2, and M phases. The equations read

dy 1 dt = ν 1b (y 7 + c) k 1b (1 + ( y 3 k 1i ) p ) + y 7 + c -k 1d y 1 , (2.2.1)
dy 2 dt = k 2b y q 1 -k 2d y 2 -k 2t y 2 + k 3t y 3 , (2.2.2)
dy 3 dt = k 2t y 2 -k 3t y 3 -k 3d y 3 , (2.2.3) dy 4 dt = ν 4b y r 3 k r 4b + y r 3 -k 4d y 4 , (2.2.4) dy 5 dt = k 5b y 4 -k 5d y 5 -k 5t y 5 + k 6t y 6 , (2.2.5)
dy 6 dt = k 5t y 5 -k 6t y 6 -k 6d y 6 + k 7a y 7 -k 6a y 6 , ( 2
∂ t n i (t, x) + ∂ x n i (t, x) + K i (y, z)n i (t, x) = 0, (2.2.11
)

n i+1 (t, 0) = ∞ 0 K i (y, z)n i (t, x)dx, ( 2
.2.12)

for i = 1, 2, and

n 1 (t, 0) = 2 ∞ 0 K 3 (y, z)n 3 (t, x)dx, (2.2.13)
The variable n i (x, t) represents the density of cells in phase i (phase 1: G1; phase 2: S/G2; phase 3: M). The variable x represents the time spent by a cell in a phase. The parameter K i is the transition rate from phase i to the next phase. The transition between phase i = 3 and phase i = 1 marks the cell division, which accounts for the coefficient 2 in the boundary condition for n 1 . Each transition rate K i depends on an average molecular state of the cells. The molecular state is given by the coupled systems of ODEs for the circadian clock and the cell cycle (Equations 2.2.1-2.2.10). The functional form of the transition rates is a Goldbeter-Koshland function [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF]:

K(y, z) = 2yJ i z -y + zJ a + yJ i + (z -y + zJ a + yJ i ) 2 -4yJ i (z -y) . ( 2

.2.14)

This function has been used to generate a switch behavior [START_REF] Novak | Mathematical model of the cell division cycle of fission yeast[END_REF]. If the ratio y/z becomes larger than one, the function switches to the upper state and the transition occurs. J a and J i are two constants that determines the stiffness of the switch, if they tend to zero, the switch tends to a step function. The transition rate from G1 (i = 1) to S/G2 (i = 2) is switched ON when the concentration of MPF reaches a certain threshold value θ 1 that instructs the cell to start DNA synthesis (K 1 = K(z 8 , θ 1 )). The G2 to M (i = 3) transition rate depends on the balance between MPF and WEE1. The cell is blocked in S/G2 and cannot transit to mitosis until MPF concentration exceeds that of WEE1 (K 2 = K(z 8 , z 9 )). The transition from M to G1, and cell division, occur when the activity of MPF goes back to baseline level

(K 3 = K(θ 3 , z 8 )
). The total cell number in each phase is given by

N i (t) = ∞ 0 n i (t, x)dx, ( 2

.2.15)

i = 1, 2, 3 and the total cell number is

N (t) = 3 i=1 N i (t).

Results

Entrainment properties

We first studied the influence of the coupling strength between the circadian clock and the cell cycle. The coupling describes the BMAL1/CLOCK-mediated rate of WEE1 activation (parameter C in equation 9). Cell cycle durations reported in the literature range from around 8 h for fast dividing lymphocytes to more than 60 h for slow tumor cells [107]. The characteristic division times of most mammalian cells coincide with the 24 h period of the day. To see how cells could entrain to the circadian clock period, we chose a cell cycle with an autonomous period (period without coupling to the circadian clock) close, but not equal to 24 h. The cell cycle period was set by scaling the time in the cell cycle equations to obtain the right period. This means that all kinetic events (activation and deactivation) are scaled uniformly. We simulated the influence of the coupling on a cell cycle with an autonomous period of 18 h, for different coupling strengths C = 0, 0.5, 1, 1.5. We observed that when the coupling strength increases, the period of the cell cycle increases (Figure 2.2). This behavior was expected for two reasons: (i) WEE1 blocks the cell cycle in G2 phase, and hence slows it down, and (ii) the period of entrainment of the circadian clock is longer than the autonomous period of the cell cycle.

We then asked whether coupling to the circadian clock always slows down the cell cycle, or whether it could speed it up. To answer this question, we looked at the influence of the coupling strength on the cell cycle for autonomous cell cycle period ranging from 8 to 60 h. Our simulations led to different modes of locking between the circadian clock and the cell cycle. For certain combinations of coupling strength and autonomous periods, the cell cycle can entrain to the circadian clock with a rational period ratio, referred to as n:m phaselocking or entrainment (Figure 2.3A). For a n:m locking, the cell divides n times each m days. These regions of the coupling strength/autonomous periods are called Arnold tongues [START_REF] Pikovsky | Synchronization[END_REF]. Arnold tongues show that the cell cycle can phase-lock to a wide range of orders with (A) Each region corresponds to an order of entrainment. The order n:m means that cells divide n times every m days. Hence, the 1:1 phase-lock region contains all cell cycles that are entrained to one division per day, or a cycle of 24 h, under an appropriate coupling strength. (B) The 24h/T' vs T plot for a fixed coupling strength (C = 1.2) has a characteristic shape, the devil's staircase [START_REF] Pikovsky | Synchronization[END_REF]. T' is the period after entrainment by the circadian clock and T is the autonomous period of the cell cycle. m up to 5. 1:1 and 1:2 phase-locks have the widest range of entrainment, but other ratios can be found for large coupling strengths, such as 2:1, 3:2 and 2:3.

For a fixed coupling strength (C = 1.2), the graph of the domains of entrainment leads to a devil's staircase (Figure 2.3B). The devil's staircase shows the frequencies (in number of cell cycles per day) of the phase-locked cell cycles as a function of the autonomous period. Phase-locked frequencies are distributed below (Figure 2.3B, shaded region) and above the autonomous frequencies (white region), indicating that entrainment by the circadian clock can either slow down, or speed up the cell cycle. The cell cycle is accelerated for intervals of autonomous periods above 24 and 48 h. Therefore, although in our model the circadian clock only acts as a break for cell cycle progression, cells with autonomous periods above 24 or 48 h can cycle faster under circadian entrainment.

Effects of coupling on growth rate

The simulations so far show that the circadian clock could make the cell cycle model run faster or slower, depending on its autonomous period. How does this translate into a net growth rate in a dividing cell population is unclear, and recent analytical results have shown counter-intuitive effects of periodic forcing on growth rates of proliferating cells. There is no systematic inequality when comparing growth rates of a population under circadian control versus a population with a constant, average control [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF][START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF], but it seems that populations under circadian control that have a cell cycle period close to multiples of 24 h proliferate faster [START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF].

We would expect the cell cycle period to be inversely proportional to the growth rate, as in the devil's staircase (Figure 2.3B). If this were so, the knowledge of the clock-entrained period should be enough to determine the cell population dynamics, without the need of population models. To test that hypothesis, and examine the effect of coupling on the growth rate, we made simulations with and without coupling to the circadian clock. Based on the Arnold tongues for C = 1.2, the autonomous period of the cell cycle was set to 20 h, inside the 1:1 phase-lock region (Figure 2.3A). When coupled to the circadian clock, the activity of MPF and WEE1 is well entrained and follows a rhythm of 24 h (Figure 2.4A,B). Driven by the new rhythm of MPF and WEE1, the fraction of dividing cells follows a rhythm of 24 h (Figure 2.4C). Even though the coupling slows down the cell cycle, the population growth rate stays practically unchanged (Figure 2.4D). This can be justified by the fact that not all cells divide at each cycle. Indeed, there was 0.38 cell division per cell per cycle with coupling, while there was 0.32 cell division per cell per cycle without coupling. Therefore, a longer cell cycle can be compensated by a larger number of division at each cycle, resulting in a higher growth rate than would be inferred from the cell cycle duration only.

To gain more insight on this non-intuitive result, we examined the impact of the coupling strength on the growth rate, for autonomous cell cycle periods ranging from 8 to 60 h. In absence of circadian coupling, the growth rate decreases almost linearly with the cell cycle period (Figure 2.5A, solid line). In presence of circadian coupling, the growth rate is decreased for most of the autonomous cell cycle periods (Figure 2.5A, non-solid lines). A notable exception is the interval between 20 h and 31 h, where the growth rate is elevated compared to the growth rate without coupling. This interval corresponds to the range of 1:1 phase-lock (Figure 2.3). For the larger coupling strengths (C ≥ 1.2), the growth rate is almost constant on this interval, as is to be expected from a synchronized population. The elevated growth rate in phase-locked populations is not systematic. For autonomous cell cycle periods above 31 h, which include the 1:2 phase-lock region, the growth rate is almost constant. In this phase-lock region, the cell cycle is entrained on a 48 h period, and the growth rate is close to the autonomous growth rate at 48 h. These results are in agreement with previous theoretical studies made with population models which showed that under circadian forcing, the growth rate was elevated near 24 h. Taken together, these results show that the growth rate is related to the entrainment of the molecular cell cycle, but that it is not possible to compare the growth rates with or without coupling [START_REF] Clairambault | Comparison of Perron and Floquet eigenvalues in age structured cell division models[END_REF][START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF][START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF].

To test the robustness and genericity of these results, we performed the same simulations on the effect of the coupling strength with two other published models for the circadian clock, one by Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF], and the other by Leloup et al. [START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF]. For the Leloup et al. model, we used parameter set 4. Both models showed the same qualitative result for the impact of the coupling strength on the growth rate. The coupling to the circadian clock increases the growth rate for periods around 24 h, over 48 h and decreases it elsewhere (Figure 2.5B,C).

Circadian clock and cancer

To investigate the role of the circadian clock in tumor development, we looked at the effect of mutations or deletions of circadian genes on the growth rate. Different types of mutations were examined, namely Per2, Bmal1, Cry2 mutations and Per2 /Cry2 double mutations. Per2 and Bmal1 mutations abolished circadian clock rhythmicity, while Cry2 and Per2 /Cry2 mutations maintained rhythmicity, in agreement with experimental data (Table 2.1 and reference [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF], details on simulating mutants are below). Two cases were studied, one considering an autonomous period of the cycle equal to 28 h and one equal to 20 h. By choosing an appropriate coupling strength (C = 1.2 for example) to the circadian clock, these two cycles could be entrained to a 1:1 cycle (Figure 2.3). We looked at changes that occurred after simulating a mutation in the circadian clock. Finally, to have a more global view, we investigated the effect of mutations for autonomous periods ranging from 8 to 60 h.

We studied the effect of mutating Per2 gene by considering that PER2 is a main ac-66 (A) tor in the negative feedback loop and simulated Per2 mutation by decreasing the rate of PER2/CRY complex formation (we set k 2b = 0.01). Simulations showed that Per2 mutants have a slower division cycle (Figure 2.6A,B, dash-dotted lines). This mutation tends to increase the growth rate in the case of 20 h autonomous cell cycle and keeps it almost equal to that of wild type cells for autonomous period of 28 h (Table 2.2). Even though the cell cycle becomes much slower in mutants, in the case of 28 h autonomous period, 1.05 cell divisions occur during each 56 h-cycle, compared to 0.4 divisions per cycle in wild type cells. This means that for Per2 mutants, some cells must divide more than once during the cycle and explains why the mutants proliferate at the same rate as the wild-type. We also supposed that PER2 activates Bmal1 transcription and simulated Per2 mutation by decreasing Bmal1 transcription rate (we set ν 4b = 1.5). Similarly, we obtained that Per2 mutants have a slower division cycle and an increased growth rate (Figure 2.7 dot-dashed lines). We simulated Bmal1 knockout by setting the transcription rate of Bmal1 ν 4b equal to 0. Simulations showed that this mutation tends to slow the cell division cycle for both 20 and 28 h autonomous period (Figure 2.6A,B, dotted lines). We observed that this mutation decreases the growth rate for autonomous periods of 28 h and increases it for autonomous periods of 20 h (Table 2.2).

0 0.5 1
We simulated deficient Cry2 mutants by decreasing the strength of the negative feedback loop (the constant k 1i was increased to k 1i = 0.8). Cry2 mutation preserves the periods of mitotic divisions (Figure 2.6A,B grey-dashed lines), but decreases the growth rate for both 20 and 28 h-autonomous the cell cycle periods (Table 2.2).

We simulated Per2 /Cry2 double mutants by assuming that Per2 mutation decreases the rate of PER2/CRY complex formation, and that Cry2 mutation decreases the strength of the negative feedback loop (we set k 1i = 0.8 and k 2b = 0.01). Our simulations showed that these double mutants have recovered a mitotic division cycle similar to that of wild type. They also have the same growth rate (Figure 2.6A,B, dashed lines; Table 2.2).

Finally, we compared the growth rates for mutants and wild type cells for autonomous cell cycle periods ranging between 8 and 60 h. Per2 mutation generally increases the growth rate. Cry2 mutation decreases it, and Bmal1 mutation increases it for autonomous periods less than 21 h and decreases it elsewhere. Our simulations also predict that Per2 /Cry2 double mutants recover a normal proliferation rate and have approximately the same growth rate for all autonomous periods of the cell cycle (Figure 2.7).

We tested the robustness and the genericity of these results by performing the same simulations on the effect of circadian genes mutations with the models proposed by Mirsky [START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF] did not show a difference in growth rate compared to the wild type (Figure 2.8). These results may be explained by the fact that this model is relatively robust to parameter variations. For the parameter set 4 in Leloup et al. [START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF], the circadian clock was most sensitive to parameters related to Bmal1, for which the effect on the population growth rate was similar to the current model and the model by Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF]. Taken together, these results predict a differential effect of certain clock gene mutations, depending on the autonomous cell cycle period of the cell population. For instance, Per2 mutant populations grow faster when the autonomous period is shorter than 40 h, but can also grow more slowly if their autonomous period is longer. Other mutations, such as Cry2, systematically slow down the population growth rate.

These results can be explained by looking at the impact of mutations on the MPF/WEE1 dynamics, which dictates the rhythm for cells to enter into mitosis and then divide. Mutationinduced change in BMAL1/CLOCK dynamics, either its period or concentration, directly influences WEE1 activity and the cell cycle dynamics. For example, in case of a 28 h autonomous period, even though the cell cycle becomes longer for Per2 mutants, the growth rate does not change.

Per2 mutation produces an arrhythmic clock with low BMAL1/CLOCK concentration (Figure 2.9B, dash-dotted line). A comparison of MPF/WEE1 dynamics between mutants and wild type cells shows how the transition rates for the M phase differ (Figure 2.10). For Per2 mutants, the transition rate is at a high level for a longer time. This means that even though the cell cycle is longer, much more cells will have the time to enter M phase and divide. This explains why growth rates are similar in Per2 mutants, even though the cell cycle period is longer. For Cry2 mutants, growth rate decreases. Cry2 mutation leads to higher rates of BMAL1/CLOCK, which in turn increases WEE1 activity (Figure 2.9B, grey-dashed line). Increasing WEE1 activity will decrease the activity of MPF (Figure 2.9A, grey-dashed line), which means that cells are blocked in G2 phase for a longer time and are prevented from transiting into mitosis. A comparison of the transition rate between wild type and Cry2 mutants shows that transition rate for wild type cells stays on a high level for a longer time than for mutant cells (Figure 2

Discussion

New insight on the regulation of the cell cycle by the circadian clock

We developed a combined molecular/population mathematical model to study how the coupling of the circadian clock to the cell cycle, through the protein WEE1, affects a proliferating cell population. The model has the novelty of combining both, intracellular and population levels. We investigated the influence of coupling on the period of the molecular cell cycle and on the growth rate of the population. The molecular model displays wide ranges of entrainment to the circadian clock, where there is a n:m ratio in the number of cell cycles and the number of circadian oscillations. We found that molecular information about the cell cycle was not always sufficient to predict how the growth rate in a dividing cell population is affected. The combined molecular/population could predict an increase in growth rate in Per2 mutants that could not be explained by the molecular model alone. We used the combined model to look at the influence of circadian clock gene mutations on the population growth rate. We found a differential effect of clock gene mutations, depending on the autonomous cell cycle period of the cell population.

Two coupled oscillators

We examined the influence of coupling the cell cycle to the circadian clock on the number of cell cycle divisions per day. We showed that for certain combinations of coupling strength and autonomous periods, the cell cycle can entrain to the circadian clock with a rational period ratio, referred to as n:m phase locking or entrainment. These regions in the coupling strength/autonomous periods space are the Arnold tongues already introduced (Figure 2.3A). Gérard and colleagues characterized domains of entrainment to 24 and 48 h periods [START_REF] Gérard | Entrainment of the mammalian cell cycle by the circadian clock: Modeling two coupled cellular rhythms[END_REF], which correspond to 1:1 and 1:2 phase-locking, respectively. While their model and the current model both predict wide ranges of 1:1 and 1:2 entrainment, there are small differences. In the current model, the 1:2 entrainment region is larger than the 1:1, while the converse is true for the model by Gérard and Goldbeter ([61], their Figure 4B, our Figure 2.3A). In the current model, 1:2 entrainment is observed at large coupling strength in a autonomous period range where 3:2 occurs (around autonomous periods of 36 h). It looks like that in the Gérard and Goldbeter model, the 1:1 entrainment takes over for larger coupling strengths. Recently, Feillet et al. used multispectral imaging of single live cells and mathematical modeling to investigate how the temporal organization of cell division at the single cell level produces daily rhythms at the population level [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF]. They demonstrated that there are multiple coexisting robust oscillatory dynamical states of the coupled clock and cell cycle in proliferating mammalian cells, namely 1:1, 5:4 and 3:2 phase locking states. We have characterized a wide range of entrainment modes, including high order phase-locking (3:5, 4:5, 5:4, 5:3), which is consistent with the experimental data of Feillet et al. and which may partly explain the observed quantized cell cycle times discussed previously [START_REF] Zamborszky | Computational analysis of mammalian cell division gated by a circadian clock: Quantized cell cycles and cell size[END_REF]. The devil's staircase provides a way to predict the frequency of cell divisions as a function of the autonomous cell cycle period. Although the circadian clock acts as a break by activating the inhibitor WEE1, for certain autonomous periods, the cell cycle frequency can still be higher with the clock than without (Figure 2.3B).

Modulation of population growth rate by the clock

We investigated the influence of coupling on the growth rate for autonomous periods of the cell cycle varying from 8 to 60 h. Clairambault and colleagues [START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF] showed, using population models, that there is no general inequality between growth rates with and without coupling to the circadian clock. Bernard et al. [START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF] found that cells under circadian control that have an interdivision time close to multiples of 24 h proliferate faster. Here, we showed that coupling increases the growth rate for autonomous periods of the cell cycle around 24 h and above 48 h. For most other periods, the growth rate is decreased. These results could not have been obtained based only on the molecular model, which predicted a smaller growth rate for autonomous periods just below 24 h.

Effect of mutating clock genes on the growth rate

We investigated the effect of single or double circadian clock gene mutations on a cell population growth rate. Fu and colleagues showed that loss of Per2 functions increased tumor development [START_REF] Fu | The circadian gene Period2 plays an important role in tumor suppression and DNA damage response in vivo[END_REF]. The roles of PER2 in the circadian clock mechanism have been unclear. It is usually considered as a main actor in the negative feedback loop, repressing the activity of BMAL1/CLOCK through the complex PER2/CRY. But some studies also suggest that PER2 activates Bmal1 transcription in an indirect manner [START_REF] Shearman | Interacting molecular loops in the mammalian circadian clock[END_REF][START_REF] Yu | Interacting feedback loops within the mammalian clock: BMAL1 is negatively autoregulated and upregulated by CRY1, CRY2, and PER2[END_REF][START_REF] Akashi | A positive role for PERIOD in mammalian circadian gene expression[END_REF]. We first examined Per2 mutation by considering that PER2 plays a repressive role in the negative feedback loop. Our simulations are in agreement with experimental results and show that Per2 mutation increases the growth rate for a wide range of autonomous periods of the cell cycle. When a positive action of PER2 on Bmal1 was assumed, we also obtained an increased growth rate for this mutation. We examined Cry2 mutation. We showed that Cry2 mutation decreases the growth rate for almost all periods of the cell cycle. This may explain the experimental results obtained by Matsuo and colleagues, who showed that the weight of regenerating liver in Cry deficient mice was significantly lower than in wild type mice [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF].

We examined Per2 /Cry2 double mutation. We showed that Per2 /Cry2 double mutants recover normal proliferation rates and have similar growth rates for all autonomous periods of the cell cycle. Oster et al. showed that inactivation of Cry2 gene in Per2 mutant mice restored circadian rhythmicity as well as normal clock gene expression patterns [119]. They showed that both the period and amplitude of Bmal1 (also of Per1 and Cry1 ) expressions in Per2 /Cry2 double mutant animals were comparable to those of wild types. Hence, if both period and amplitude of Bmal1 are comparable to those of wild type, WEE1 profile will not be changed for these double mutants, preserving normal dynamics for the cell cycle.

We also explored Cry1 mutation and Cry1 /Cry2 double mutation, based on the assumption that CRY1 plays a more important role in the negative feedback loop [START_REF] Langmesser | Interaction of circadian clock proteins PER2 and CRY with BMAL1 and CLOCK[END_REF]. In the current model, these mutations completely abolished cell proliferation (data not shown). There is no experimental evidence that disruption of the circadian clock can totally prevent cell cycle progression, and it is likely that the cell cycle relies on factors not included in the current model to proceed through division.

Robustness of the results

To test the robustness of our results, we performed our main simulations with two other models for the circadian clock (Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF], Leloup et al. [START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF]). We simulated the effect of coupling and circadian genes mutation on the growth rate. Both models show the same qualitative result that we obtain for the impact of the coupling on the growth rate. Results on mutating circadian genes obtained with the model proposed by Mirsky et The current circadian clock model [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF] shows limitations in reproducing experimental data that are inherent to models with simplifying assumptions and distinct molecular species lumped together. To check how the simplifying assumptions affect the results, we made a detailed comparison between the current model and the model by Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF] on the effects of circadian gene mutations. What it is critical in our model is the effect specific mutations on the activity of BMAL1/CLOCK (period and concentration), which regulates directly WEE1 and the cell cycle. Hence, we simulated the effects of circadian gene mutations on BMAL1/CLOCK using Mirsky et al. model (in the same way the authors did it in the original study, Table S3 in [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF]) and compared the results with those obtained by the current model. Both models showed similar effects for Bmal1 and Per2 mutations. The main differences between the two models are the effects of Cry2 and Per2 /Cry2 mutations. Mirsky et al. predict a longer period for Cry2 mutants, which may be more realistic in the case of lung explants and fibroblasts. The model by Mirsky et al. predicts an arrhythmic clock for Per2 /Cry2 double mutation, in contrast to the current model and experiments showing normal rhythmicity for these double mutants [119] . The model by Forger and Peskin also predicts a rhythmic clock for this double mutation [START_REF] Forger | A detailed predictive model of the mammalian circadian clock[END_REF]. The main limitation of the current model is the way Cry1 mutation and Cry1 /Cry2 double mutation are approached. Both are simulated in the same way, by decreasing further the strength of the negative feedback loop. This leads to an arrhythmic clock with high constitutive BMAL1/CLOCK concentration inhibiting cell proliferation, which cannot be supported by experimental data (data not shown). Results given on the growth rate by the model of Mirsky et al. on Cry1 mutation seems more coherent. Cry1 mutants have a decreased growth rate compared to wild types for a wide interval of autonomous periods of the cell cycle (Figure 2.8A, grey solid line).

Conclusion

Combining a molecular model to a population model offers new insight on the influence of the circadian clock on the growth of a cell population. Disruption of the circadian clock can increase or decrease the growth rate, as well as the period of mitotic divisions, depending on which clock gene is affected. In some cases, even though the cell cycle slows down, the growth rate can still increase, making the combination of a molecular model and population model unavoidable to study the effect of circadian clock disruption. This can have beneficial impacts on chronotherapy, which aims to develop new strategies in cancer therapies by a better understanding of the circadian clock and its impact on cell proliferation.

The combined model presented in this study is the first step in developing a fully multiscale model for the interaction between the circadian clock and the cell cycle. The multiscale model describes a cell population p structured with a molecular content (y, z) describing the circadian clock and the cell cycle. Heterogeneity among cells can be fully taken into account in a multiscale model, but at the cost of a high-dimensional phase space (here 10D). Even though it has limitations, the current molecular model is simple enough to be amenable to a multiscale description, which, in our view, is essential.

Several studies have shown that the tolerance and the toxicity of drugs varies according to their administration time [START_REF] Focan | Circadian rhythms and cancer chemotherapy[END_REF][START_REF] Hrushesky | The application of circadian chronobiology to cancer chemotherapy[END_REF][START_REF] Lévi | Therapeutic implications of circadian rhythms in cancer patients[END_REF]. Clinical studies showed that compared to standard chemotherapies, chronomodulated chemotherapies, which aim to deliver drugs at an optimal time of the day, could be more efficient and better tolerated by patients [START_REF] Focan | Continuous delivery of venous 5-fluorouracil and arterial 5-fluorodeoxyuridine for hepatic metastases from colorectal cancer: Feasibility and tolerance in a randomized phase II trial comparing flat versus chronomodulated infusion[END_REF][START_REF] Focan | Chronotherapy with 5fluorouracil, folinic acid and carboplatin for metastatic colorectal cancer; an interesting therapeutic index in a phase II trial[END_REF]. In a recent work, Bernard et al. [START_REF] Bernard | Tumor growth rate determines the timing of optimal chronomodulated treatment schedules[END_REF] used a simple cell population model under chronomodulated treatment and developed a quantitative method to identify biological parameters important for the successful design of a chronotherapy strategy. They found that optimal times depend not only on the circadian status but also on the cell cycle kinetics of the tumor. They suggested that the length of the cell cycle is important to determine the best treatment times and intervals. For fast growing tumors, with short S phase, administering a drug that targets the S phase of the cell cycle at 28.8 h intervals may be safer than treating at 24 h intervals, and that for slow growing tumors, with a long S phase, treating at 24 h intervals would be the best option. The circadian clock is often disrupted in advanced stage cancers, perhaps because this gives a competitive advantage to growing tumor cells [START_REF] Savvidis | Circadian rhythm disruption in cancer biology[END_REF]. The current model could be useful to predict how the cell cycle is modified following circadian clock disruption. Combined with the method proposed by bernard et al. [START_REF] Bernard | Tumor growth rate determines the timing of optimal chronomodulated treatment schedules[END_REF], this may be of great importance to determine the right time for drug delivery. By taking into account complex interactions between the cell cycle, the circadian clock and the treatment, the combined molecular/population model can be a helpful tool for chronotherapy.

Chapter 3 A particle method for high-dimensional transport equations: application for the regulation of the cell cycle by the circadian clock 3.1 Introduction

In the previous chapter we have presented a model that couples the circadian clock and the cell cycle through transition coefficients in an age-structured equation. Though this model had the novelty of capturing the influence of intracellular dynamics on the growth rate of cells, it lacks a multiscale description and cannot take into consideration, for example, intracellular heterogeneity among cells. In this chapter we present a multiscale version of the previous model and a numerical method for solving it. This is done by structuring our transport equation by the molecular contents of the coupled cell cycle-circadian clock oscillator. This means that instead of taking a transport equation of the form

∂ t ρ(a, t) + ∂ a ρ(a, t) = f (x)ρ(a, t),
where a is the time spent by cells in a phase of the cycle and f (x) is some function that depends on the intracellular components x = (x 1 , • • • , x d ) as in chapter 2, we now take an equation of the form

∂ t ρ(x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = L[x, λ] ρ(x, t, λ) .
Partial derivatives for the convective term in this latter equation are now taken with respect to the molecular contents x and the space where we solve this equation is a d-dimensional space with d the number of molecular components x i ; this is why we call it a molecularstructured equation. Our model becomes high-dimensional, with d typically larger than 10, and classical numerical methods such as finite volumes/differences are not appropriate for solving this transport equation. We circumvent this difficulty by using a particle method.

Particle methods have arisen as an alternative to classical numerical methods for solving high-dimensional problems. They are used in different applications, for example the incompressible Euler equations in fluid mechanics [START_REF] Harlow | The Particle-in-Cell Method for Fluid Dynamics[END_REF][START_REF] Leonard | Vortex methods for flow simulation[END_REF][START_REF] Leonard | Computing three-dimensional incompressible flows with vortex elements[END_REF], the Vlasov equation in plasma physics [START_REF] Bossy | A stochastic particle method for the McKean-Vlasov and the Burgers equation[END_REF][START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF] and in turbulence models for reactive flows [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF][START_REF] Pope | PDF methods for turbulent reactive flows[END_REF]. The computational implementation of particle methods is conceptually simple. At a given time, the solution is represented by a large number of particles, each having its own properties, for example position and weight. These properties evolve in time according to a system of stochastic or deterministic differential equations and the classical solution to the PDE is obtained by reconstructing the particles distribution. Theoretically, the solution is thought as a linear combination of Dirac masses

ρ(x, t) = N j=1 α j (t)δ(x -x j )
where α j is the weight of particle j, x j its position (state) and N is the total number of particles. To obtain a classical solution, one has to update the positions and weights of particles and then regularize the Dirac masses. The overriding strength of particle methods is that, for N fixed, the size of the system increases only linearly with the dimensionality of the space. This means that if we have a structured equation with dimension d, we have to solve d × N ODEs/SDEs. Since recovering a classical solution needs a regularization of the dirac distribution, the performance of particle methods depends also on the quality of the regularization procedure. For generalities about particle methods, one can refer to [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF][START_REF] Cottet | Vortex methods: theory and practice[END_REF].

Solving our equation using the particle method becomes natural, since each particle can be associated to a single cell. We make several improvements on the model presented in Chapter 2. Here, we consider a cell population that is heterogeneous, which means that each cell has its own properties, namely different molecular concentrations and intrinsic cell cycle period. We assume that cells are coupled to the same master circadian clock and that cells are connected together by assuming that the rate of Per/Cry mRNA production of each cell is dependent on the average Per/Cry production of all cells. As for the previous model, we consider that division is dependent on the MPF-WEE1 antagonistic activities. In addition, we assume that cell growth is dependent on the total number of cells. Specifically, we consider a limited growth by assuming that MPF activity decreases with increasing number of cells. These additional assumptions, and the fact that the transport equation is structured by molecular concentrations confer to the model its multiscale nature.

This chapter is divided into three parts. In the first one, we introduce the structure of our model and give details about the equations we use. In the second part, we introduce deterministic particle methods and their mathematical background. In the third part, we expose a test case for the particle method and the results of our model. In addition we have developed a code (in C) that is adaptable to different intra-cellular dynamics. This means that the code we present in the appendix can be used to solve similar hyperbolic type equations with different convective velocity u.

Description of the model

Equation without cell birth or death

We have a large collection of cells with state (x, λ) The distinction between the dynamical and the parameter state is that the dynamical state changes during the lifespan of a cell, while the parameter state is fixed over the lifespan of the cell, but can vary among cells or when a cell divides. The population is described by its density ρ(x, t, λ). The evolution of density is described by a nonlinear, nonlocal first-order hyperbolic equation

∈ Ω × Λ ⊂ R d × R p ,
∂ρ ∂t (x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = 0. (3.2.1)
We assume that the dynamics of the cellular state x is governed by a deterministic system of ODEs, represented by the term u(x, t, λ, ψ), that depends on the cellular state x and λ and on m population-level statistics ψ : Ω × R × Λ → R m where

ψ i = ρ, Φ i = Ω×Λ ρ(x, t, λ)Φ i (x, t, λ)dxdλ, i = 1, • • • , m (3.2.2) with Φ i : Ω × R × Λ → R, i = 1, • • • , m taken such that ψ i is finite.

Equation with cell birth and death and parameter variation

Equation 3.2.1 with a source term has the following general form

∂ρ ∂t (x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = L[x, λ] ρ(x, t, λ) . (3.2.3) L[x, λ
] is an operator describing the population relative growth rate. We impose that L be only determined by the cellular state (x, λ) and not by the population statistics ψ. The fact that the growth rate does not depend on the population statistics ψ is voluntary, as cells rely on their internal state only to determine their fate.

We seek solution of the form ρ(x, t, λ) =

N (t) j=1 α j (t)δ (x j (t),λ j (t)) .
The approximate solution ρ lies in a measure space M(R d × [0, T ] × R p ). When computing the solution, we can choose to aggregate several cells in single states, or keep all cells at their individual states, or use a mixture of aggregation and individual-state description. We favored the individual-state description. This implies that α j ≡ 1 for all i = 1, • • • , N, and that the total cell number is N (t). The solution has the form

ρ(x, t, λ) = N (t) j=1 δ (x j (t),λ j (t))
.

We distinguish between two kinds of growth rates, deterministic and stochastic. For the deterministic growth rate, cell divides when (x, λ) crosses an open subset of a d -1 dimensional smooth manifold Γ ∈ Ω, while for the stochastic growth rate, cells divide with a finite rate r(x, λ). For the cellular state to cross Γ, the characteristic lines must be transversal to Γ, i.e.

u(x, t, λ) • ∇F (x) = 0,

where F defines implicitly Γ, with F (x) = 0 for all x ∈ Γ. In particular u(x, t, λ) must not vanish on Γ, and the crossing must always occurs in the same direction: u(x, t, λ) • ∇F (x) has the same sign in all connected subsets of Γ. Γ lies in Ω but can depend on λ as parameter (Figure 3.1). When the cell crosses Γ the growth rate is given by

r(x, λ) = δ (x∈Γ) .
In practice, defining Γ that satisfies the transversality condition may not be practical, so we can add a condition that the crossing should be in a specific direction: u • ∇F > 0 or u • ∇F < 0. Cell death is treated in a similar way. In the deterministic case, because there is a finite number of cells N (t), divisions occur at discrete times. At division time t, N (t) = N (t -) + 1, and the parameter state of the daughter cells are updated, according to a probability density p(λ|z) of jumping from state z to state λ.

For the stochastic growth rate, the relative growth rate r(x, λ) is given by a finite function of (x, λ), with the interpretation that the probability for a single cell to divide in a small time interval t is t × r(x, λ). The parameter state can be thought of as a differentiation marker, and we assume that changes in λ occur only at cell division. The dynamics on λ is based on [START_REF] Bernard | How to build a multiscale model in biology[END_REF]. Equation 10 in [START_REF] Bernard | How to build a multiscale model in biology[END_REF] describes the evolution of a cell density in a differentiation state space, analogous to the parameter state space here. The evolution equation is an integral-differential equation of the form

∂ρ ∂t (t, λ) = Λ p(λ|z)R(z)ρ(t, z)dz -R(λ)ρ(t, λ) + r(λ)ρ(t, λ). (3.2.4)
The kernel p(λ|z) is the probability density of jumps from parameter state z to state λ, R is the rate at which jumps are occurring (the differentiation rate) and r is the growth rate of the population. Because change in parameter λ occurs only at division, the differentiation rate is linked to the to growth rate. For instance, if there is no cell death, then R = r, i.e. differentiation occurs only when there is a cell division. The full equations on (x, λ) now reads

∂ρ ∂t (x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, tλ) = Λ p(x, λ|z)R(x, z)ρ(x, t, z)dz -R(x, λ)ρ(x, t, λ) + r(x, λ)ρ(x, t, λ). (3.2.5)
When the jumps in the parameter space are not to large and the dependence of p on z not too important, we can approximate the integral term by a diffusion. Let σ 2 (x, z) be the variance of the jump distribution at z. Then the equation can be approximated as

∂ρ ∂t (x, t, λ) + ∇ x • u(x, t, λ, ψ)ρ(x, t, λ) = 1 2 Δ λ • σ 2 (x, t, λ)R(x, λ)ρ(x, t, λ) + r(x, λ)ρ(x, t, λ). (3.2.6)

Intracellular dynamics

Equation (3.2.1) describes the evolution of cells in the space of intracellular components. We are interested in studying the effect of circadian regulation on the cell cycle. To do that, we take the intracellular dynamics to be the coupled circadian/cell cycle deterministic ODE system that we presented in the previous chapter

dx dt = u(x, t, λ, ψ). (3.2.7)
where

u 1 (x) = dx 1 dt = ν 1b (x 7 + c) k 1b (1 + ( x 3 k 1i ) p ) + x 7 + c -k 1d x 1 + k s ψ 1 , (3.2.8) u 2 (x) = dx 2 dt = k 2b x q 1 -k 2d x 2 -k 2t x 2 + k 3t x 3 , ( 3
.2.9)

u 3 (x) = dx 3 dt = k 2t x 2 -k 3t x 3 -k 3d x 3 , (3.2.10) u 4 (x) = dx 4 dt = ν 4b x r 3 k r 4b + x r 3 -k 4d x 4 , (3.2.11) u 5 (x) = dx 5 dt = k 5b x 4 -k 5d x 5 -k 5t x 5 + k 6t x 6 , ( 3
.2.12) We recall that equations (3.2.8-3.2.14) describe the circadian oscillator and (2.2.8-2.2.10) the cell cycle. Dynamical variables are x 1 , Per2 or Cry mRNA and proteins; x 2 , PER2/CRY complex (cytoplasm); x 3 , PER2/CRY complex (nucleus); x 4 , Bmal1 mRNA; x 5 , BMAL1 cytoplasmic protein; x 6 , BMAL1 nuclear protein; x 7 , Active BMAL1; x 8 , Active MPF; x 9 , Active WEE1; x 10 , Active MPF inhibitor. Coupling of these two oscillators is taken into consideration in equation (3.2.16) through the regulation of WEE1 by BMAL1/CLOCK (in the term Cx 7 , C is the coupling strength). Details about this model can be found in the previous chapter. The cellular dynamics is now part of a multiscale system, some changes were added to the system to reflect that cells interact with each other. These changes are detailed in the coming paragraphs. Peripheral tissues are not known to couple their circadian clocks. This example of coupling in proliferating cells is taken as an illustration of the multiscale modeling.

u 6 (x) = dx 6 dt = k 5t x 5 -k 6t x 6 -k 6d x 6 + k 7a x 7 -k 6a x 6 , (3.2.13) u 7 (x) = dx 7 dt = k 6a x 6 -k 7a x 7 -k 7d x 7 , (3.2.14) u 8 (x) = dx 8 dt = λ k lmpf +k 0mpf exp -ηψ 2 k n 1mpf k n 1mpf +x n 8 +sx n 10 (1 -x 8 ) -d wee1 x 9 x 8 , (3.2.15) u 9 (x) = dx 9 dt = λ k actw k actw + d w1 (c w + Cx 7 ) + k actw k actw + d w1 -1 k inactw x n 8 x 9 k n 1wee1 + x n 8 -d w2 x 9 , ( 3

Population synchronization

Including connectivity between cells is in certain tissues important to maintain a robust synchronized activity. It is known, for example, that neuronal clocks within the the suprachiasmatic nucleus SCN form a heterogenous network that must synchronize to maintain time keeping activity. The coherent output of the SCN is established by intracellular signaling factors, such as vasointestinal polypeptide [START_REF] Enright | Temporal precision in circadian systems: A reliable neuronal clock from unreliable components?[END_REF]. A simple way to induce synchronization in our model is to make Per/Cry mRNA transcription depends on the average concentration of Per/Cry among cells. For that, we add to equation (3.2.8) the term k s ψ 1 := k s ρ, Φ 1 with k s a coefficient that represents the connectivity strength and Φ 1 = x 1 . Peripheral tissues are not known to couple their circadian clocks. This example of coupling in proliferating cells is taken as an illustration of the multiscale modeling.

Variability among cells

Variability among cells arises naturally from the difference in their molecular contents. This is taken into consideration in our model through the initial states of cells, which can be chosen randomly. We add another source of variability through the parameter λ. We assume that each cell has a distinct cell cycle period. This is modeled by multiplying equations (3.2.15-3.2.17) by a scaling factor λ. Cell division can also induce variability by assigning to each new born cell an intrinsic cell cycle period that is different from its mother cell. This is taken into consideration by taking a non trivial distribution p(λ|z).

Cell division

An important aspect of our model is that it takes into account cell division. Cell division is based on the antagonistic relationship between the mitosis promoting factor MPF and the protein WEE1. It is assumed that a cell enters mitosis once the activity of MPF surpasses that of WEE1 and then divides once MPF activity shuts down abruptly. Based on this mechanism we consider two types of division, one deterministic and one stochastic. The deterministic division occurs "exactly" every time MPF activity rises above WEE1 activity and then shuts down. This was taken into consideration by considering a growth rate r(x, λ)

:= r d (x, λ) = δ (x∈Γ) .
For the stochastic division, we consider that a cell divides with a certain probability t × r(x, λ) within a small time step t. The division is a function of x and λ that mimics the deterministic case. For example, a switch-like function that takes small values on one side of Γ and large values on the other side could be used. Here we used the Koshland-Goldbeter switch defined in the previous chapter (see equation 2.2.14).

Limited growth

To make sure the growth is bounded (for physiological and computational reasons), we assume that cell proliferation stops when the total number of cells reaches a threshold value. We assume that the activation coefficient of MPF is decreasing at an exponential rate proportional to the total number of cells. This is taken into consideration in equation (3.2.15) by multiplying MPF activation coefficient k 0mpf by exp -ηψ 2 where ψ 2 = ρ, Φ 2 with Φ 2 = 1 (ψ 2 = N (t)). When cell number is large enough, MPF activity cannot increase above that of WEE1, and cell division is blocked.

Description of the particle method

The theoretical results we present in this section are based on the work of P.A. Raviart [START_REF] Raviart | An analysis of particle methods[END_REF].

A d-dimensional population equation, in its conservative form, can be written as follows

∂ρ(x, t) ∂t + d i=1 ∂ ∂x i (u i (x, t)ρ(x, t)) + u 0 (x, t)ρ(x, t) = f (x, t), x ∈ R d , t > 0 (3.3.18) ρ(x, 0) = ρ 0 (x). ( 3
φ ∈ C 1 0 (R d × [0, T [) (space of continuously differentiable functions with compact support in R d × [0, T [) T 0 R d ρL * φdxdt = T 0 R d fφdxdt + R d ρ 0 φ(•, 0)dx
where L * is the adjoint of the differential operator L defined by:

Lρ = ∂ρ ∂t + d i=1 ∂ ∂x i (u i ρ) + u 0 ρ and L φ = - ∂φ ∂t - d i=1 u i ∂φ ∂x i + u 0 φ.
The particle method consists of considering even less regular solutions, hence finding a solution of equation (3.3.18) in the space of measures of subsets on R d . By definition, a

function ρ ∈ M(R d × [0, T [) is a measure solution of (3.3.18) if ρ, L φ = f, φ + ρ 0 , φ(•, 0) ∀φ ∈ C 1 0 R d × [0, T [ (3.3.20)
where M(S), for a set S, is the space of measures defined on S and •, • is the duality pairing between M(S) and C 0 0 (S). Equality (3.3.20)

makes sense if ρ 0 ∈ M(R n ) and f ∈ M(R n × [0, T [).
The intuition behind the particle method is to start with a distribution of particles that approximate the initial condition and then follow the evolution in time of the positions and weights of these particle according to the velocity u and the functions u 0 and f . Since the particle solution is a measure solution, this implies it is irregular. To obtain a solution in the usual classical sense at a give time T , one has to recover the classical solution with regularization techniques. Hereafter, we describe the different steps of obtaining the solution for equation (3.3.18) using the particle method. For the sake of simplicity we consider f = 0.

Step 1: approximation of the initial condition

Given an initial condition of our problem ρ 0 ∈ C 0 (R d ), we would like to chose the initial distribution of our particles x j and their weights α j so that ρ 0 h = j α j δ(x-x j ) approximates ρ 0 . This should be understood as an approximation in the sense of measures, which means that we are looking for a test function φ ∈ C 0 0 (R d ) such that

ρ 0 , φ = R d ρ 0 φdx and ρ 0 h , φ = j α j φ(x j ).
This yields a typical numerical quadrature problem where one would like to approximate R d ρ 0 φ by j α j φ(x j ). A simple, but not necessarily most efficient, way for choosing α j and

x j can be obtained by using the midpoint quadrature rule. If we partition our domain R d into a collection of sets B j , and place the particles at t = 0 into the centers of mass of each set: the second-order midpoint quadrature rule leads to the following initial set of particle distribution:

x j = x j c , α j = μ(B j )ρ 0 (x j c
), where x j c and μ(B j ) are respectively the center of mass and the area of set j. Supposing our domain Ω is a d-cube, and dividing it into regular sets of size h such that x j = jh, the approximation of the initial density reads

ρ 0 h (x) = j h d ρ 0 (x j )δ(x -x j ). (3.3.21)
This approximation means that if ρ 0 ∈ C 0 (R d ), then we have that for all functions

φ ∈ C 0 0 (R d ) lim h→0 ρ 0 h -ρ 0 , φ = 0. (3.3.22)
The proof of this convergence result in M(R d ) can be deduced from the following theorem concerning quadrature formulas in R d . Theorem 5. [START_REF] Raviart | An analysis of particle methods[END_REF] Let m ≥ 1 be an integer and let p > d m , q = p p-1 , then there exists a constant C > 0 independent of h such that for all function

g ∈ W m,p (R d ) ∩ L 1 (R d ) if m ≤ 2 or g ∈ W m,p (R d ) ∩ W m-1,1 (R d ) if m ≥ 3, we have R d g(x)dx -h d j g(x j ) ≤ Ch m+ d q j | g | m,p,B j
where W m,p (R d ) is the Sobolev space of functions which are, together with their deriva-

tives up to order m in L p (R d ) and | g | m,p,Ω := |α|=m ∂ α g p L p (Ω) 1/p is the usual Sobolev semi-norm.
Using 5, we get that lim

h→0 | R d g(x)dx -h d j g(x j ) |= 0 for any function g ∈ C 0 0 (R d ).
Finally the convergence result 3.3.22 follows by taking g = ρ 0 φ.

Step 2: updating positions and weights

The particles positions X j (t) (to not confuse between X j and x j , our notations implies that X j (0) = x j ) and their corresponding weights α j (t) can be computed at a given time t by solving the following ordinary differential system:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ d dt X j (t) = u(X j (t), t) d dt α j (t) = -u 0 (X j (t), t)α j (t)
X j (0) = x j and α j (0) = α j and the measure solution at a given time t, is given by

ρ h (x, t) = j α j (t)δ(x -X j (t)). (3.3.23)
Step 3: classical solution recovery Solution (3.3.23) is a set of particles with the positions in Ω and weights. To obtain a solution in the classical sense, one has to regularize the particle solution. This can be done by taking the convolution product of ρ h with a cut-off function

ζ ∈ C 0 (R d ) ∩ L 1 (R d ) parameterized with such that R d ζ(x)dx = 1 and ζ (x) = 1 d ζ( x )
. The classical solution is hence defined as

ρ h (•, t) = ρ h (•, t) * ζ = j α j (t)ζ (x -X j (t)). (3.3.24)
Convergence of the regularized particle solution to the classical one is obtained using the following L p estimate for the error ρ(•, t)ρ h (•, t) with 1 ≤ p ≤ +∞ Theorem 6. Assume the following conditions hold

• there exists an integer k ≥ 1 such that (i) R d ζ(x) = 1 (3.3.25) (ii) R d x α ζ(x)dx = 0 ∀α ∈ R d and 1 ≤ α ≤ k -1 (3.3.26) (iii) R d | x | k | ζ(x) | dx < +∞ (3.3.27) • the cut-off function ζ belongs to the space W m,∞ (R d ) ∩ W m,1 (R d ) for some integer m > d.
• the coefficients

u i ∈ C 0 (R d × [0, T ]) and u 1 , • • • , u d , u 0 + div u ∈ L ∞ (0, T, W l,∞ (R d ))
where l = max(k, m).

Then if ρ 0 ∈ W l,p (R d ), there exists a constant C = C(T ) > 0 such that for all t ∈ [0, T ] ρ(•, t) -ρ h (•, t) L p (R d ) ≤ C k ρ 0 k,p,R d + ( h ) m ρ 0 m,p,R d ) . (3.3.28)
We give here a brief description of the proof of the above theorem. For more details one could refer to [START_REF] Raviart | An analysis of particle methods[END_REF]. The idea consists of splitting the error ρ(•, t)ρ h into a sum of two terms u(•, t)u(•, t) * ζ + ρ(•, t)ρ(•, t) * ζ and giving an L p bound for each one. We start by stating the following lemma Lemma 7. Assuming there exists a constant k such that conditions 3.3.25-3.3.27 are verified, then we have for some constant C > 0 and for all function

f ∈ W k,p (R d ), 1 ≤ p ≤ +∞ f * ζ -f L p (R d ) ≤ C k | f | k,p,R d Now if ρ 0 ∈ W k,p (R d
), it follows from the smoothness of the coefficients u i that ρ belongs to L ∞ (0, T ; W k,p (R d )). Hence applying lemma 7 we get an L p bound for the first term of the sum

ρ(•, t) -ρ(•, t) * ζ L p (R d ) ≤ c k | ρ(•, t) | k,p,R d ≤ c k ρ 0 k,p,R d . (3.3.29)
The constant c needs not to be the same in both inequalities, however to simplify the notations we kept it the same.

To show the L p bound for the second term of the sum we remark that

((ρ(•, t) -ρ h (•, t)) * ζ )(x) = R d ρ(y, t)ζ (x -y)dx - j α j (t)ζ (x -X j (t)) = j E j (g(x, •, t))
where E j (g) = B j g(x)dxh d g(x j ) and g(x, y, t) = ρ 0 (y) exp(-t 0 u 0 (X(s; y, 0)ds)ζ (x -X(t; y, 0)).

After checking that y → g(x, y, t) ∈ W 1,m (R d ), we use theorem 5 with p = 1 to obtain

| ((ρ(•, t) -ρ h (•, t)) * ζ )(x) |≤ ch m | g(x, •, t) | m,1,R d .
Now using the smoothness of the coefficients u i we get

| g(x, •, t) | m,1,R d ≤ ch m |α|+|β|≤m R d | ∂ α ρ 0 (y)∂ β ζ (x -X(t; y, 0)) | dy.
Remarking that the right hand side of the above inequality is a finite sum of a the convolution products

| ∂ α ρ 0 (X(t; •, 0)) | * | ζ | and using the fact that f * g L r (R d ) ≤ f L p (R d ) + g L q (R d ) with 1 r = 1 p -1 q -1, we get that | g(x, •, t) | m,1,R d ≤ ch m |α|+|β|≤m ∂ α ρ 0 (X(t; •, 0)) L p (R d ) ∂ β ζ L 1 (R d )
and by remarking that

∂ β ζ L 1 (R d ) ≤ c |β| and that ∂ α ρ 0 (X(t; •, 0)) L p (R d ) ≤ ∂ α ρ 0 L p (R d )
we finally get

((ρ(•, t) -ρ h (•, t)) * ζ )(x) ≤ c h m ρ 0 m,p,R d . (3.3.30)
Therefore the proof of the needed L p bound 3.3.28 is obtained by adding the L p bounds inequalities 3.3.29 and 3.3.30 which proves theorem 6.

Results

2-d test problem

To test the particle method and our code, we consider an advection problem that was studied in [START_REF] Durran | Numerical Methods for Wave Equations in Geophysical Fluid Dynamics[END_REF][START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF]. In this example, a passive tracer is advected in a non divergent deformational flow. The spatial domain is the unit square [0, 1] × [0, 1] and the initial condition is a cosine bell. The equations read ∂ρ(x, y, t) ∂t

+ ∂ ∂x (u 1 (x, y, t)ρ(x, y, t)) + ∂ ∂y (u 2 (x, y, t)ρ(x, y, t)) = 0 (3.4.31) ρ 0 (x, y) = 1 + cos(πr(x, y)) 2 , ( 3.4.32)
where r(x, y) = min 1, 4 (x -1/4) 2 + (y -1/4) 2 , and the velocity field is a swirling shear flow defined as it follows:

u 1 (x, y, t) = sin 2 (πx) sin(2πy) cos(πt/5) (3.4.33) u 2 (x, y, t) = -sin 2 (πy) sin(2πx) cos(πt/5), (3.4.34)
The tracer distribution is most deformed at t = 2.5. Since the velocity periodically reverses direction, the tracer distribution comes back to its initial position after a period T = 5. We compared the solution obtained with the particle method and with a classical nonoscillatory method (Third order-upwind for the spatial derivative and third order stability input : Number of cells N Index of WEE1 and MPF i W and i M Division type: Dirac or Koshland or no division time step τ output: A matrix S that contains the state of cells at a given time T t = 0 while t <= T do t = t + τ % Advance one time step using 4 th order Runge-Kutta scheme:

for j = 1 to N do for i = 1 to M do k (1) i = u(t, x j i , λ j , ψ) k (2) i = u(t + τ 2 , x j i + k 1 2 τ, λ j , ψ) k (3) i = u(t + τ 2 , x j i + k 2 2 τ, λ j , ψ) k (4) i = u(t + τ, x j i + k 3 τ, λ j , ψ) S i,j = S i,j + τ 6 (k (1) i + 2k
(2) i + 2k

(3) i + k (4) i ) end end for j = 1 to N do flag j 1 = 1 (S i M ,j >S i W ,j ) flag j 2 = 1 (S i M ,j <S i W ,j ) × 1 (S i M ,j <θ 1 ) K G = G(S i W ,j , S i M ,j ) generate a random number r uniformly distributed in [0, 1] if division type = Dirac then if (flag j 1 = 1 and flag j 2 = 1) then N = N + 1, λ N +1 = λ j for i = 1 to M do S i,N +1 = S i,j , end flag j 1 = flag j 2 = 0 end end else if division type = Koshland then if (flag j 1 = 1 and r < K G × τ ) then N = N + 1, λ N +1 = λ j for i = 1 to M do S i,N +1 = S i,j , end flag j 1 = flag j 2 = 0 end end else
Go back to the main time loop end end end Algorithm 1: Simplified Pseudo algorithm to solve equation 3.2.5 using the particle method. The full code is given as an Appendix. In this algorithm the molecular contents and the intrinsic cell cycle period are kept the same after division. This can be changed by assigning for example random values to S i,N +1 and λ N +1 after division.

preserving Runge-Kutta for time derivative [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]) at both t = 2.5 and t = 5. To make a reasonable comparison, we took 10000 particles for the particle method and a uniform 250 × 250 grid for the classical one. The solution using the particle method is irregular compared to the classical solution (Figure 3.2A,B). This is to be expected due to the nature of particle methods which consider solutions in the measure sense. Specific regularization techniques can be applied to regularize the solution such as the filtering algorithms for oscillatory solutions of hyperbolic PDEs presented in [START_REF] Engquist | Nonlinear filters for efficient shock computation[END_REF]. Another possibility is to take larger values of for the cut-off function ζ . However, increasing the value of makes the solution more diffusive (Figure 3.2C,D). This can be overcome by modifying the convolution procedure, and taking a cut-function ζ i with variable width that depends on the distance between particles [START_REF] Chertock | On a practical implementation of particle methods[END_REF]. Our interest in this study is to examine statistic-like quantities such as molecular averages and synchronization properties. Such information can be deduced from the particles distribution without the need of a "best regular" solution recovered from them as it is the case in other applications like fluid dynamics. Therefore, we do not focus on advanced reconstruction techniques for the solution, because the particles distribution, even though not regular, allows us to obtain information we are looking for. Also, the fact that the particle method enables us to reduce our 10-d problem to a system of ODEs remains the unavoidable advantage of this method. Another advantage of particle methods over classical non-oscillatory schemes is their shape preserving property over time. This is clearly illustrated by comparing the solutions at t = 5 where we see how the classical scheme leads to some numerical diffusion in the solution while the particle solution remains accurate (Figure 3.2E,F). This is important in our study, since we are following proliferation of cells, hence for long time intervals, we do not lose information about cells state.

Regulation of the cell cycle by the circadian clock

We solve equation (3.2.3) with L = 0, which means without division. We assume that each cell has an intrinsic cell cycle period included in the interval [20 h, 28 h] . This is done by assigning to each particle a different coefficient λ. We let all cells have the same initial molecular state and consider two cases; one without coupling to the circadian clock (C = 0), and one with coupling to the circadian clock (C = 1.2). We plot a projection of the solution in the subspaces (x 7 , x 9 ) and (x 7 , x 8 ), which are the subspaces of molecular concentrations (BMAL1/CLOCK, WEE1) and (BMAL1/CLOCK, MPF). We obtain a random distribution of particles for the non-coupled case and a limit cycle distribution for the coupled case (Figure

3.3)

. This means that the circadian clock is forcing all cells to oscillate with the same period. This is in agreement with the results obtained by the model presented in the previous chapter.

We have shown that a population of cells with cell cycle periods between 20 h and 28 h are brought to oscillate with a unique cell cycle period of 24 h for large values of C.

Cell division

We took initially 4 cells, two of them with λ = 230 and two with λ = 180 (these values are taken randomly and do not necessarily lead to a cell cycle period of integer value). We solve equation (3.2.3) with a source term and consider two cases, one with r = r d and the other with r = r s . For both cases we first take η = 10 -16 , which means that MPF activity does not depend on the total cells number (taking a very small value for η insures that the coefficient exp -ηψ 2 in equation (3.2.15) is equal to one, implying no dependence of k 0mpf on the total cell number). We assume that new born cells retain the same molecular concentrations of their mother cells at the division time, hence the differentiation rate R is equal to 0. This insures that MPF cycle does not change along birth. Our simulations show that the total number of cells after four days is equal to 1152 for the first case, 900 for the second one with a stiff Koshland-Goldbeter function and 480 with a non stiff one. This is natural and is justified by the fact that the division rate r d is a Dirac-type division which means that division depends deterministically on MPF activity. Looking at MPF activity, we see that it peaks 9 times for λ = 230, and 6 times for λ = 180 every 4 days (Figure 3.4A,B). If division occurs exactly once MPF activity accomplish a normal cycle, the total number of cells should be 2 × 2 9 + 2 × 2 6 = 1152, which is the result obtained with r = r d (Figure 3.4C). The division rate r s introduces stochasticity in the decision for division. In this case, division does not depend only on MPF-WEE1 activity, but also on the probability r s × t at each time step t Figure (3.4C). Second, we assumed that proliferation depends on the total number of cells, and took η = 2 × 10 -3 (increasing the value of η will make the factor exp -ηψ 2 less than one and hence decreases the value of MPF activation coefficient k 0mpf ). We took initially 100 cells with different intrinsic cell cycle periods and followed the total number of cells over 15 days. We obtain that after 8 days, division stops and the curve reaches a steady state (3.4D). Increasing the value of η decreases the activity of MPF with increasing number of cells, which means that at a certain time, MPF activity will decrease below a threshold that does not allow the cell to proceed through mitosis.

Synchronization

To study the synchronization of cells, we compared the molecular concentrations of a cell chosen randomly with the average molecular concentrations for all cells. We let the rate of production of Per2/Cry mRNA of each cell to be dependent on the average value for Per2/Cry mRNA. This was done by taking k s = 0.05 in equation (2.2.1). We took a population of 100 cells with autonomous cell cycle periods randomly distributed between 20 and 28 hours. Initial molecular concentrations are chosen randomly between 0 and 1 and each cell has a distinct initial molecular state. We did not consider division, hence L = 0 and we do consider coupling between the cell cycle and the circadian clock (C = 1.2). We followed the evolution of 3 components, Per2/Cry mRNA, BMAL1/CLOCK and WEE1, over 20 days. Our simulations show that for k s = 0, the average molecular concentrations tend to have small oscillations which are asynchronous with the oscillations of the random cell concentrations (Figure 3.5A,B,C). Whereas for k s = 0.05, we obtain that the average molecular concentrations tend to coincide with the concentrations of the random cell (Figure 3.5E,F,G). This indicates that all cells are oscillating in a similar manner (same period and phase), and indicates synchronization of all cells of the population. We would like to emphasize here that, in the case k s = 0, even though cells are coupled in the same manner to the circadian clock, cells keep oscillating in an asynchronous matter. The circadian clock regulate all cells to oscillate in a similar period which is 24 h in this case, but due to the randomness in the initial molecular concentrations, each cell oscillate with a different phase.

Heterogeneity of cells and growth rate

We studied the growth rate of a cell population where each cell has a different cell cycle period. We took initially 100 cells and confer to each of them a coefficient λ chosen randomly between 128 and 193. These coefficients scale the intracellular cell cycle system so that periods range randomly from 12 h to 18 h. The circadian control strength value C was fixed to 1.6, we did not consider connection between cells and we did not consider dependence of the intracellular dynamics on any population level statistics. Simulations were done over 20 days and with a net death rate to mainly reduce the computation time. Our simulations showed a growth rate with two peaks per day, indicating that there are two rounds of division every day (Figure 3.6A). To explain this bimodal behavior for the growth rate, we looked at the average activity of WEE1 and MPF. Simulations show that MPF activity overcomes that of WEE1 once a day (Figure 3.6B), suggesting that division should occur only once a day. This cannot explain the two peaks per day obtained for the growth rate. This underlies the fact that biological function markers that most of the time relies on average-like information can mislead the interpretation in some cases. For that, we investigated the individual-like division mechanism. We followed a subpopulation of cells and examined at which time each cell divided. There are almost 3 regimes for division: 1 division per day, 2 divisions per day and 3 divisions every 2 days (Figure 3.6C,D). This is not surprising and can be explained by the entrainment results obtained in the previous chapter for autonomous cell cycle periods between 12 and 18 h. The double peaks for the growth rate can be then explained by the fact that all cells are dividing at least once a day and some of them two times per day or

Discussion and conclusion

We presented in this chapter a multiscale model for the regulation of the cell cycle by the circadian clock. The model is structured by the molecular contents of the coupled circadiancell cycle oscillator and is high dimensional. We presented a solution using a deterministic particle method and developed a code that is adaptable to different molecular systems. Several studies, similar to ours, dealt with modeling cell population connected to molecular systems influencing cell growth. Bekkal Brikci et al. developed a nonlinear model for the dynamics of a cell population divided into a proliferative and quiescent compartments. This model was structured by the time spent by a cell in a proliferative phase and by the amount of Cyclin D CDK4/6 [START_REF] Brikci | An age-and-cyclin-structured cell population model for healthy and tumoral tissues[END_REF][START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF]. Ribba et al. developed a multiscale model of cancer growth based on the genetic and molecular features of the evolution of colorectal cancer. This model investigated the role of gene-dependent cell cycle regulation in the response of tumors to therapeutic irradiation protocols [START_REF] Ribba | A multiscale mathematical model of cancer, and its use in analysing irradiation therapies[END_REF]. In more recent works, Prokopiou et al. presented a multiscale computational model to study the maturation of CD8 T-cells in a lymph node controlled by their molecular profile [START_REF] Prokopiou | Multiscale modeling of the early CD8 T-cell immune response in lymph nodes: An integrative study[END_REF]. Even though these models contain details about both, molecular and population level, they are of low dimensionality. Our model captures more information but has the advantage of being of dimension 10. As a consequence, we had to find a way to simulate transport equations in high-dimensions at reasonable cost. The particle method does just that. In addition the code we developed is adaptable to any system of differential equations for the advective term. This means that the code can give solutions for even higher dimensional equations. Also, none of these models has taken into consideration heterogeneity among cells; which is an important feature of our model. There exist such type of models, aiming at studying the synchronization of a population of cells. Namely, the work by Kunz and Acherman who studied, using van der Pol oscillators, inter cellular coupling mechanisms between generic oscillators and showed how locally coupled networks can robustly synchronize [START_REF] Kunz | Simulation of circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators[END_REF]. Rougemont et al. used more abstract Kuramoto oscillators, in which only the phase is described, with periods and phases randomly varying in time to characterize the source of phase dispersion [START_REF] Rougemont | Collective synchronization in populations of globally coupled phase oscillators with drifting frequencies[END_REF]. As examples of models using realistic genetic networks to describe synchronization, we cite the work by Roenneberg and Merrow who proposed the concept of Zeitnehmer, where the cellular circadian oscillator feeds back on the input pathway of Zeitgebers, blurring the distinction between intra and A extracellular components [START_REF] Roenneber | Demasking biological oscillators: Properties and principles of entrainment exemplified by the neurospora circadian clock[END_REF]. A more comprehensive study was developed by Bernard et al. who proposed a heterogenous network of circadian neuronal oscillators model for the synchronization of circadian oscillators that combines intercellular and single cell-level dynamics [START_REF] Bernard | Synchronization-induced rhythmicity of circadian oscillators in the suprachiasmatic nucleus[END_REF]. This category of models can inform about coupling properties of a cell population but cannot take into consideration cell proliferation.

Globally our model captures information on cell proliferation, intra-and intercellular dynamics. The use of particle methods in our study was of crucial importance. It allowed us to take into consideration a large system describing molecular dynamics and structure the transport equation by many variables. In such contexts, the method is natural and simple since each cell can be associated to a particle. This can pave the way for the use of high-dimensional structured models since particle methods, as we have shown, makes the resolution of high dimensional equations possible.

Chapter 4 Conclusion

In this thesis, we developed a multiscale mathematical model for the regulation of the cell cycle by the circadian clock. We studied the influence of a circadian control on the proliferation of a cell population.

We started our study by a simplified version of the multiscale model. We developed a model that consists of two coupled systems, one for the proliferation of cells and one for the molecular cell cycle/circadian clock connected network. Cell proliferation was described by a system of age-structured transport equations with transition coefficients that depend on the molecular state of cells. We studied the way the cell cycle entrains to the circadian clock and characterized multiple domains of entrainment. We showed that there are two wide ranges of entrainment 1:1 and 1:2, and other ranges which are narrower, for example 2:1, 3:2, 5:4 and 3:4. Such regimes of entrainment were found in recent experimental studies, for example in the work of Feillet et al. [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF] who obtained 1:1, 3:2 and 5:4 entrainment modes. Sandler et al. showed that cells follow a 2:1 entrainment to the circadian clock [START_REF] Sandler | Lineage correlations of single cell division time as a probe of cell-cycle dynamics[END_REF]. We showed that the circadian clock can slow down or accelerate the cell cycle depending on cell autonomous cell cycle period. We studied the influence of a circadian control on the net growth rate of cells. We showed that the circadian clock increases the growth rate for autonomous cell cycle periods around 24 h and larger than 48 h. We examined the impact of clock gene mutations on the growth rate and showed that disruption of the circadian mechanism leads to abnormal cell proliferation. Namely, we showed that Per2 mutation increases the growth rate, Cry mutations decreases it, and Bmal1 mutation can either decrease it or increase it depending on the autonomous cell cycle period. This is in agreement with experimental studies, namely the study of Fu el al. [START_REF] Fu | The circadian gene Period2 plays an important role in tumor suppression and DNA damage response in vivo[END_REF] who showed that Per2 have increasing sensitivity to tumor development and Matsuo et al. who obtained impaired liver regeneration in Cry deficient mice [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF]. This model has the novelty of combining both, a population and a molecular description. It gave new insights on the regulation of the cell cycle by the circadian clock that could not be obtained by a population model or a molecular model alone.

Next, we presented the multiscale model. It consists of a transport partial differential equation structured by the molecular concentrations of the coupled cell cycle-circadian clock system. The main difficulty encountered with this model is its high-dimensionality aspect. We used a deterministic particle method to circumvent this difficulty and developed a code (in C) to solve the partial differential equation. This code can be adapted to different convective velocities (intracellular dynamics). We endowed this model with several additional properties that confers to it its multiscale nature; namely, the heterogeneity among cells and the dependence of the intracellular dynamics on population like statistical quantities, such as the total number of cells or a total average concentration. We examined the solution of the structured equation considering different cases. We showed that without birth and with a circadian control, the particle solution follows a well defined trajectory in the state space, unlike without circadian control where particles stay randomly distributed. This yields a different view for entrainment properties than the one studied using the simplified model. We examined a case where the rate of MPF production depends inversely on the total cell number and showed how this stops cell proliferation. We showed how the dependence of the rate of Per/Cry production on the total average Per/Cry concentration for all cells induces synchronization properties. We finally examined the growth rate of a cell population with cells having different autonomous cell cycle periods. We obtained two peaks of division per day that could not be explained by average molecular concentrations and explained it by following divisions time for each cell. This highlights the importance of our modeling approach, since it captures both, population dynamics and individual-like mechanisms.

Our model can be helpful for cancer chronotherapy. Chronotherapy aims to take benefit of circadian rhythms to find an optimal time for drug delivery. This makes the drug less toxic and more efficient. Several studies showed that the cell cycle dynamics determines optimal treatment time. We have seen how a clock gene mutation leads to abnormal proliferation and different dynamical behaviors for the fraction of cells in mitosis. Disruption of circadian rhythms is often associated with advanced stage tumors. Our model can help estimate the time cancer cells are dividing and hence helps finding the appropriate time for the drug delivery. Also, two peaks growth rate, like the one we obtained with our multiscale model, were observed in mice tumors [START_REF] Wood | Circadian clock coordinates cancer cell cycle progression, thymidylate synthase, and 5-fluorouracil therapeutic index[END_REF]. It was shown that the circadian clock component BMAL1 and protein WEE1 of tumor cells vary throughout the day in synchrony with tumor growth, mitotic index, 5-Fu drug target and the toxic therapeutic index of 5-FU. Such results can be reproduced by our model and hence make it a useful tool for chronotherapy. The work presented in chapter 4 is under preparation for submission.

Several improvements can be added to our modeling as a future continuation of this work. Improvements can be done on both, the intracellular and population levels. On the intracellular level, we considered for the circadian clock the model developed by Becker-Weimann et al. More detailed models can be considered such as the ones developed by Mirsky et al. [START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF] or Forger et al. [START_REF] Forger | A detailed predictive model of the mammalian circadian clock[END_REF]. Our description for the cell cycle was limited to the antagonistic activity of MPF and WEE1. A more comprehensive model for the cell cycle can also be considered. This makes it possible to take into consideration more than one link between the cell cycle and the circadian clock. Also, we note that we considered only a one-way coupling, which means an implication for the circadian clock on the cell cycle. However, several studies showed that also cell cycle components regulates circadian ones. This also can be further considered. The multiscale description yields a high dimensional model; this constrained our choice for a reduced system on the molecular level. However, as we indicated in chapter 4, the code we developed can be adapted to different convective terms u. Hence, a future work can be done by considering a more elaborate system that describes the coupling between the cell cycle and the circadian clock. On the population level, we considered that cell division depends only on the state variables, more types of division can be considered. The particle method and the code should be improved to be able to take into consideration a more general source term L. Improvements on the multiscale description can be inspected, namely on the mechanisms inducing cells heterogeneity and connectivity. Our description for heterogeneity is limited to variability among molecular concentrations and cell cycle periods. Its is known, for example, that cells circadian mechanism is not the same for all cells in some tissues. Connection through different signaling pathways is found to insure circadian synchronization among these tissues. Such kind of variability and connectivity can be added to improve our modeling.

During the thesis, we have also worked on other projects that are not directly related to the circadian clock/cell cycle coupling modeling.

We developed a mathematical model (presented in Appendix A) to estimate meningioma tumor age for several patients who had a surgery to extract the tumor. In this work we analyzed, using minimization techniques, experimental data about cells average age, Ki-67 index, net growth of cells related to patients age, and tumor volume. We were able to fit the data and give an estimate for the time the tumor started growing for each patient. We are examining additional data about tumor volume, and once the analysis is done, this work will be the subject of a publication.

We also worked on a project about modeling chronic myelogenous leukemia treatment with the Tyrosine Kinase Inhibitor, Imatinib. Briefly, in this work we had clinical data for the BCR-ABL ratios for a group of 104 patients who were monitored during Imatinib therapy. We developed a mathematical model based on a system of ordinary differential equations that describe the proliferation of a population of leukemic and autologous immune cells. We observed that patients leukemic load was able to partially or fully suppress the autologous immune response. We, therefore, defined an immune window(a range of leukemic loads) for which the autologous immune system induces an improved response. We deduced that Imatinib therapy drives the leukemic load under this immune window allowing patients immune cells to expand and eventually mount an efficient recognition of the residual leukemic burden. Our study led us to hypothesize that immunotherapy may complement Tyrosine Kinase Inhibitors treatments, by helping to maintain a patients autologous immune response when the leukemia stimulus alone is insufficient. This work was concretized by a publication that was recently submitted (BCR-ABL transcripts variations in chronic phase chronic myelogenous leukemia pa- tients on imatinib: Possible role of the autologous immune system. Geoffrey D. Clapp, Thomas Lepoutre, Raouf El Cheikh, Samuel Bernard, Jérémy Ruby, Hélène Labussière-Wallet, Franck E. Nicolini, Doron Levy).

Circadian rhythms disruption is linked to several diseases. Understanding their implication on our health is important. We focused on understanding the implication of circadian rhythms on the proliferation of a cell population. We developed a multiscale model for this purpose. Our approach is novel and innovative. While most of the previous works in the same theme limited their studies to one aspect, such as molecular dynamics or population dynamics; our multiscale model describes different aspects, namely cell proliferation, individual-like dynamics and cells connectivity.

Describing complex biological processes needs advanced mathematical modeling tools. Multiscale modeling is an unavoidable choice. This thesis makes a step forward in this direction.

Globally the thesis work was mostly focused on the regulation of the cell cycle by the circadian clock (two published papers and one under preparation) and two side projects; one on modeling chronic myelogenous leukemia (one submitted paper) and one on modeling meningioma tumor growth (one paper under preparation). patient: the tumor volume and the average age of tumor cells. A third quantity, the growth coefficient, was estimated from literature [START_REF] Nakamura | The natural history of incidental meningiomas[END_REF]. Using the model and available data, we estimated the values of three parameters for each patient: a proliferation coefficient describing how fast the tumor grows initially, a carrying capacity describing the maximal size the tumor can reach, and the time elapsed since the beginning of tumor development (the time tumor had a volume corresponding to one cell). Parameters were estimated using a least-square method that minimizes the difference between collected data and model predictions. Our simulations lead to non intuitive results for the relation between the tumor size and its age.

A.1.2 Model description

Gompertz model describes the evolution of the number of cells in time with the following ordinary differential equation:

⎧ ⎨ ⎩ dN (t) dt = α ln K N (t) N (t), N (0) = N 0 . (A.1.1)
Here, N (t) is the number of cells at time t (years), α is the proliferation coefficient, K is the carrying capacity or the maximum amount of cells that can be attained, and N 0 is the initial number of cells. The solution of this equation can be given explicitly by

N (t) = Ke ln N 0 K e -αt .
This model describes the evolution of the tumor size, but does not give any information about the age of cells themselves. To take the age of cells into account, it is necessary to endow the Gompertz model with an age-structure so that we can investigate both, the tumor size and the age of cells at a given time. Adding an age structure x, our model reads:

∂n(x, t) ∂t + ∂n(x, t) ∂x = 0, (A.1.2a) n(x, 0) = N 0 δ(x), (A.1.2b) n(0, t) = αN (t) ln K N (t) . (A.1.2c)
Here, n(x, t) is the density of cells of age x at time t. N 0 is the initial number of cells, the Dirac delta function δ(x) insures that all cells, initially, have an age x = 0 years. The The solution of the age-structured Gompertz equation can be determined explicitly:

n(x, t) = ⎧ ⎪ ⎨ ⎪ ⎩ N 0 δ(x -t) + αK ln K N 0 e -α(t-x) e ln N 0 K e -α(t-x) if t ≥ x, 0 i f x > t, (A.1.3)
We made several assumptions regarding cell dynamics during aging. First, we assume that there is no cell death (0 on the right hand side of equation A.1.2a). Second, we assume that a cell is born at age 0, (boundary condition at x = 0 in equation A.1.2c). Third, mother cells retain their age after division. These assumptions have an impact on the age distribution of cells that will be discussed later.

A.1.3 Tumor age estimate: lack of data

The analytical solution A.1.3 gives the density of cells for a known couple (age x, time t) and given parameters (α, K). However, our goal is to find the values of (t, α, K) that give an average age for cells and a volume that are equal to those collected experimentally.

To solve this inverse problem, we define the following cost functional

J 1 t, α, K; Ṽ , ã =| Ṽ -V | 2 + | ã -a | 2 ,
where, Ṽ and ã are respectively the volume and the average age predicted by the model, V and a are respectively the volume and the average age given experimentally. Variables are scaled to reduce the optimal search interval. We took one volume unit as 10 -9 cm 3 and the initial population is set to one cell, that is N 0 = 1.

To find the tumor age, we should solve the following minimization problem. Find: min (t,α,k) J 1 t, α, K; Ṽ , ã such that (t, α, K) belong to an admissible set of values. This implies that the growth time t should not exceed the patient age, the coefficients α and K should be chosen so that the value of the growth rate α ln(K/N ) should not exceed the growth rate induced by the Ki-67 index. This is due to the fact that we do only consider net growth, and this may not reflect the heterogeneity of the population, because there might be another subset of cells which are not proliferating. Therefore, the Ki-67 index, gives us an upper bound for the growth rate.

The minimization problem is solved numerically using the Matlab function fminsearch. Simulations show that there exists an interval of values for (t, α, K) that minimize the function J 1 (Figure A.1). This infer that there is not a unique parameter set (t, α, K) that fits experimental data about the volume and the average age; this was to be expected because, roughly speaking, we have three parameters to estimate and two variables to fit.

A.1.4 Additional data

The main obstacle encountered in minimizing the function J 1 is the lack of information regarding the proliferation coefficient of tumor cells. We could not predict a unique value for the tumor age. We obtained an interval of possible values that all fit the experimental data.

To circumvent this obstacle, we relied on experimental data coming from the work by Nakamura et al. [START_REF] Nakamura | The natural history of incidental meningiomas[END_REF]. This work investigated the relation between the patient age and the meningioma growth rate. These experimental data helped us to estimate the required growth coefficient for tumor cells. Nakamura data can be helpful in two different ways: estimating the absolute growth and estimating the relative growth. We analyze both ways in the coming paragraphs.

First approach: estimating the absolute growth

In a first attempt, we looked at the absolute growth per year given by Nakamura data. It helped us get global estimates for α and K in the Gompertz model. It was possible to choose, Fitting the absolute growth data using Gompertz model. Solid curve is obtained by finding the appropriate couple values (α, K) yielding an absolute growth (V (t +1)-V (t)) that fit best nakamura data. V (t) is the volume at time t and is given by V c × N (t), where V c is the volume of one cell.

for each patient, a unique triplet (t, α, K) that is consistent with Nakamura data and obtain a fitting curve that estimates the absolute tumor growth per year (Figure A.2). Relying on this, we went back to our initial optimization problem, which yield an interval of values for (α, K), and looked at which of these values could give a curve of absolute growth per year similar to the one estimated by Nakamura data (Figure A.3). Once (α, K) were fixed, we were able to give a unique estimate for the tumor age (Table A

Second approach: estimating the relative growth

Rather than looking at the absolute growth per year, we looked at the relative growth, which is known to depend on the tumor size. We follow the same approach and fit the relative growth rate as a function of tumor size, using Gompertz model (Figure A.5). The fitted curve gives us a relative growth coefficient α ln( K N ) that depends on the tumor size or the number of cells. Assuming that, for a tumor size V (number of cells = N f ), the fitting curve yields a relative growth equal to γ r , we impose a new constraint to the cost function

J 2 =| Ṽ -V d | 2 + | ã -a d | 2 + | γr -γ r | 2 ,
where γr = α ln( K N f ). This additional constraint ensures that the cost function returns a (A)

Final approach

We took benefit of Nakamura data in two different ways to estimate the tumor age. Both led to reasonable results. We end up this investigation by combining both approaches. We define the following cost functional

J 3 =| Ṽ -V | 2 + |< ã > -< a >| 2 +δ 1 | γr -γ r | 2 +δ 2 | γa -γ a | 2 ,
where γr and γa (γ r = α ln( K N f ), γa = V (t +1)-V (t)) are respectively the relative growth and the absolute growth rates to be found. γ r and γ a are the corresponding growth rates given by Nakamura's data. δ 1 and δ 2 are weight constants that play the role of finding an average value between γ r and γ a .

Minimizing J 3 enabled us to estimate the tumor age for all patients while respecting experimental data (Table A.4).

A.1.5 Discussion

In this work, we propose a mathematical model to estimate the age of meningioma for several patients who were operated to extract the tumor. The model consists of an age-structured transport equation coupled with a Gompertz type equation to model the growth of cells. Analysis was based on three main experimental data for each patient: the tumor volume, the average age, and the growth coefficient given by the experimental work of Nakamura et al. We point out that we had an additional experimental data about the percentage of proliferating cells (Ki-67 index) for each sample. Unfortunately, the Ki-67 index did not help to find the proliferative ability of tumor cells. Even though, it informs us on the percentage of cycling cells, we were not able to use this percentage to find the growth coefficient. None of the estimated values for (α, K) could achieve a relative growth α ln( K N ) equal to the one given by the Ki-67 index. This can be explained by the fact, that the population of cells may be heterogenous, or there is a subset of cells that is not proliferating. Hence, we concluded that the growth rate induced by the Ki-67 index is an upper bound for the growth coefficient. By defining a minimization function that takes in consideration the volume, the average age and the growth coefficient (both relative and absolute), we were able to estimate the tumor age. Results are counter-intuitive in certain cases. One would intuitively expect that larger tumors are older. However, the link between the the tumor size and its age is more intricate than that. The tumor volume does not affect much its age. For instance, looking at figure A.7, we observe that patient 1 has a tumor of volume 50.5 cm3 but the estimated age is small. This can be explained by the fact that tumor cells have a big proliferation coefficient and the tumor size is rapidly increasing. Conversely, patient 6 has a small tumor (2.25 cm3) but the tumor is very old, and this is is due to the low proliferation coefficient, so tumor size is increasing slowly. Based on parameter estimates for each patient, we speculate whether the tumor is dangerous or not. If the curve is in its steady state (for example patient 1 in figure A.7), the tumor reached a limit size value and will not grow much more than that. Whereas, if the curve is still in its exponential phase, tumor is still growing and could be dangerous, --------------------------------------------- --------------------------------------------------------------------h ----------------------------------------------------------------------- * -------------------------/ ----------------------------------------------------------------------r -------------------------------------------------------------------------- * = params [START_REF] Achermann | Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves[END_REF] ; double k1b

// In c a s e you c h o o s e k o s h l a n d : S p e c i f i y t h e s t i f n e s s o f t h e s w i t c h // s m a l l e r v a l u e s g i v e s t i f f e r s w i t c h

double const K G f u n c t i o n s t i f f = 0 . 1 ; double const K G d i v i s i o n c o n s t a n t = 1 0 . ; //

// I n i t i a l s t a t e : //-----------------------------------------------// 1--In c a s e you c h o s e a f i x e d i n i t i a l s t a t e //-

//-----------------------------------------------// 2--In c a s e you c h o o s e random i n i t i a l s t a t e //-----------------------------------------------// Choose t h e range " each " i n i t i a l v a l u e t o b e l o n g " randomly " t o : // y i j b e l o n g s t o [ a , b ] ; a and b s h o u l d NOT be i d e n t i c a l and a

a u p t . h

//--------------------------// Newing f r e e i n g s p a c e //-

i g h t h a n d s i d e s a b i n a . cpp

= params [START_REF] Akashi | A positive role for PERIOD in mammalian circadian gene expression[END_REF] ; double k 1 i = params [START_REF] Altinok | An automaton model for the cell cycle[END_REF] ; double c = p a r a m s [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF] ; double p = params [START_REF] Aschoff | Tagesperiodik bei Mäusestämmen unter konstanten Umgebungsbedingungen[END_REF] ; double k1d

= params [START_REF] Aschoff | Circadian control of body temperature[END_REF] ; double k2b

= params [START_REF] Aschoff | Angeborene 24-Stunden-Periodik beim Käcken[END_REF] ; double q = p a r a m s [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF] ; double k2d

= params [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF] ; double k2t

= ; dydx [START_REF] Achermann | Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves[END_REF] [ j ] = k2b * pow ( y1 , q ) -( k2d + k2t ) * y2 + k 3 t * y3 ; dydx [START_REF] Akashi | A positive role for PERIOD in mammalian circadian gene expression[END_REF] [ j ] = k 2 t * y2 -k3t * y3 -k3d * y3 ; dydx [START_REF] Altinok | An automaton model for the cell cycle[END_REF] [ j ] = ( nu4b * pow ( y3 , r ) ) / ( pow ( k4b , r ) + pow ( y3 , r ) ) -k4d * y4 ; dydx [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF] [ j ] = k5b * y4 -(k5d + k 5 t ) * y5 + k 6 t * y6 ; dydx [START_REF] Aschoff | Tagesperiodik bei Mäusestämmen unter konstanten Umgebungsbedingungen[END_REF] [ j ] = k 5 t * y5 -(k6t + k6d + k6a ) * y6 + k7a * y7 ; dydx [START_REF] Aschoff | Circadian control of body temperature[END_REF] [ j ] = k6a * y6 -k7a * y7 -k7d * y7 ; dydx [START_REF] Aschoff | Angeborene 24-Stunden-Periodik beim Käcken[END_REF] [ j ] = s c a l f a c t o r 2 * ( ( 4 . + k0mpf * exp(-lambda * n s t a t e s ) ) * (1.0 -y8 ) * pow ( k1mpf , nmw) / ( pow ( k1mpf , nmw) + pow(y8 ,nmw) + s * pow ( y10 , nmw) ) -dwee1 * y9 * y8 ) ; dydx [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF] [ j ] = s c a l f a c t o r 2 * ( kactw / ( kactw + dw1 ) * ( cw+kbm * ( y7 -0.9629) + 0 . 9 6 2 9 ) + ( kactw /( kactw + dw1 ) -1 . ) * kinactw * pow ( y8 , nmw) / ( pow ( k1wee1 , nmw)+pow ( y8 , nmw) ) * y9 -dw2 * y9 ) ; dydx [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF] [ j ] = s c a l f a c t o r 2 * k a c t * ( y8 -y10 ) ; } MVfunction out . c l e a r ( ) ; } -------------------------------------------------------------------------- *

/ * ------------------------------------------------------------------------O D E s y s t e m p a r a m s a s s i g n a t i o n . cpp

/------------------------------------------------------------------------------// p a r a m e t e r s f o r t h e c o u p l e d o s c i l l a t o r //------------------------------------------------------------------------------

params [ 0 ] = . 0 4 2 ; // s c a l f a c t 1 params [START_REF] Achermann | Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves[END_REF] = . / params [ 0 ] ; // nu1b params [START_REF] Akashi | A positive role for PERIOD in mammalian circadian gene expression[END_REF] = . ; // k 1 b params [START_REF] Altinok | An automaton model for the cell cycle[END_REF] = . 5 6 ; // k 1 i params [START_REF] Aschoff | Zeitgeber der tierischen Tagesperiodik[END_REF] = . 0 1 ; // c params [START_REF] Aschoff | Tagesperiodik bei Mäusestämmen unter konstanten Umgebungsbedingungen[END_REF] = . ; // p params [START_REF] Aschoff | Circadian control of body temperature[END_REF] = . 1 2 / params [ 0 ] ; // k1d params [START_REF] Aschoff | Angeborene 24-Stunden-Periodik beim Käcken[END_REF] = . 3 / params [ 0 ] ; // k 2 b params [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF] = . ; // q params [START_REF] Becker-Weimann | Modeling feedback loops of the mammalian circadian oscillator[END_REF] = . 0 5 / params [ 0 ] ; // k2d params [ 1 0 ] = . ---------------------------------------------------------------------- -------------------------------------------------------------------------- * ------------------------------------------------------------------ ---------------------------------------------------------------------- * n n s t a t e s = n s t a t e s ; double * * K; // S t a t e s M a t r i x double * * dydx ; // double * * s c a l f a c t 2 ; // e q u i v a l e n t t o lambda i n t h e ODE s y s t e m i n t * * d i v c o u n t e r 1 ;

// e q u i v a l e n t t o f l a g 1 i n t h e pseudo-a l g o r i t h m i n t * * d i v c o u n t e r 2 ;

// e q u i v a l e n t t o f l a g 2 i n t h e pseudo-a l g o r i // P a r a m e t e r s i n i t i a t i o n O D E s y s t e m p a r a m s a s s i g n a t i o n (&params [ 0 ] ) ; p r i n t p o i n t e r (&params [ 0 ] , nparams ) ; f o r ( i n t i = 0; i < 1 ; i ++) { f o r ( i n t j = 0; j < n s t a t e s ; j ++) { d i v c o u n t e r 1 [ 0 ] [ j ] = i n d 1 ; d i v c o u n t e r 2 [ 0 ] [ j ] = i n d 2 ; } } // i n i t i a l s t a t e i n i t i a l s t a t e (&K[ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , nstates , i n i t i a l v a l u e s t y p e , i n i t s t a t e v a l u e s f i x e d , i n i t s t a t e v a l u e s r a n d , s c a l f a c t 2 f i x , c e l l c y c l e p e r i o d m i n f , cell c y c l e p e r i o d m a x f ) ; // P r i n t i n g i n t i t i a l pop o f s t r e a m w r i t e i n i t p o p ( " i n i t p o p . t e x t " ) ; p r i n t m a t r i x ( w r i t e i n i t p o p , &K[ 0 ] , neqs , n s t a t e s ) ; w r i t e i n i t p o p . ----------------------------------------------------------------------- * -------------------------------------------------------------------s ----------------------------------------------------------------------- * f o r ( i n t j = n s t a t e s ; j < n n s t a t e s ; j ++) {

// c o u t << Time << s e t w ( 1 2 ) << n n s t a t e s << e n d l ; // U p d a t i n g t h e s y s t e m one t i m e s t e p u s i n g

t a t e s u p d a t e . cpp

//----------------------------------------------------//----------------------------------------------------// U p d a t i n g t h e p o p u l a t i o n a f t e r each t i m e s t e p // G e n e r a t i n g new c e l l i f a c e l l p a s s e s from M t o G1 // The new c e l l has t h e same m o l e c u l a r c o n c e n t r a t i o n o f i t s mother ? // R e a l l o c a t i o n o f t h e m a t r i x

// R e a s s i g n i n g f o r e a c h c e l l t h e same p e r i o d o f // i t s mother ! can be o t h e r w i s e ! !

i f ( strcmp ( i n i t i a l v a l u e s , "random" ) == 0) { f o r ( i n t i = 0; i < n e q s ; i ++) { K

Résumé :

This thesis is dedicated to the development of a multiscale mathematical model that describes the regulation of the cell cycle by the circadian clock. What motivated this work is the fact that several tumorigenic diseases are linked to circadian rhythms disruption. We would like to understand the effect of circadian rhythms on the proliferation of a cell population and hence give plausible explanation for diseases that arise form circadian clock disruption. The mammalian cell cycle and the circadian clock are two molecular processes that operate in a rhythmic manner and exquisite precision. On one hand, the cell cycle is driven by the rhythmic activity of cyclin-dependent kinases which dictate the time a cell must engage mitosis and the time it must divide giving birth to two daughter cells. On the other hand, the circadian clock is a system of transcriptional and translational feedback-loops that generates sustained oscillations of different mRNAs and proteins with a period of approximately 24 h. It turns out that several components of the circadian clock regulates various cyclin-dependent kinases at different stages of the cell cycle. This makes the circadian clock a key player of the temporal organization of the cell cycle and makes these two biological processes act as two tightly coupled oscillators. Our modeling approach consists of using a molecular-structured partial differential equation that describes the proliferation of a cell population. Proliferation depends on the coupled cell cycle-circadian clock molecular state of cells. Due to the large number of molecular components involved in the cell cycle-circadian clock system, the problem becomes of high-dimensionality and specific numerical techniques are needed to solve the equation. As a first step, we simplify the problem, and use a system of transport partial differential equations structured by the time spent by cells in a phase of the cell cycle. This system is coupled to the molecular one, via transition coefficients that depend on the molecular state of cells. Even tough it is a simplified version, this model has the novelty of combining both population and intracellular levels. We use it to study the entrainment of the cell cycle by the circadian clock and the effects of regulation on the net growth of cells. Afterwards, we pass to the fully multi-scale model and use the particle method to circumvent the highdimensionality aspect. This method consists of representing the population of cells by a large number of particles, each having its own set of properties, position and weight. These properties evolve in time according to a system of ordinary differential equations, so that the particles simulate individual cells evolving in the molecular state space. The solution of the system can be reconstructed from individual particles and shown to converge to the solution of the initial system. The main advantage of particle methods is their usefulness for high dimensional problems where classical numerical methods as finite difference/volumes/elements fail. Key Words : circadian clock, cell cycle, dynamical systems, structured transport partial differential equations, entrainment, growth rate, particle method.

R

  . El Cheikh, T. Lepoutre and S. Bernard. Modeling Biological Rhythms in Cell Populations. Mathematical Modelling of Natural Phenomena / Volume 7 / Issue 06 / 2012, pp 107 -125.

1. 1

 1 Van der Pol oscillator: simulations with different values for the damping coefficient and bifurcation diagram . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Goodwin oscillator: simulations with different values for the Hill coefficient and Bifurcation diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Simulations showing the oscillatory activity for the components of the circadian clock with Becker-Weimann et al. and Goldbeter models . . . . . . . . . 1.4 One way and toggle switches . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Simulations showing the oscillatory activity and the limit cycle for the Albert Goldbeter cell cycle model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Molecular network scheme for fission yeast cell cycle . . . . . . . . . . . . . . 1.7 Simulations for the cell cycle model of Tyson and Novak . . . . . . . . . . . 1.8 Solution of the age-structured model with and without periodic transition coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Molecular schemes for the circadian network and for the coupling between the cell cycle and the circadian clock . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Simulations for the van der Pol socillator with different values of the damping coefficient : (A) Positive values for lead to a limit cycle. (B) Negative values for lead to damped oscillations. (C) Bifurcation diagram for the van der Pol oscillator. Solid black line: stable steady state; dashed black line: unstable steady state; blue and red lines: limit cycle.

Figure 1 . 3 :

 13 Figure 1.3: (A) Becker-Weimann et al. circadian model: proteins and mRNAs concentrations. Solid line: total BMAL1 complex (variables: y 5 + y 6 + y 7 ); dashed line: PER2/CRY protein (variable y 3 ); dot-dashed line: Per2/Cry mRNA (variable y 1 ); dotted line: Bmal1 mRNA (variable y 4 ). BMAL1 protein oscillates antiphasic to Per2/Cry mRNA and with a period approximately equal to 24 hours. PER2/CRY protein oscillates with a phase delay of 7.5 hours compared to Per2/Cry mRNA. (B) Goldbeter circadian clock model. Solid line: Per mRNA; dot-dashed line: unphosphorylated PER form; dotted line: monophosphorylated PER form; grey line: bi-phosphorylated PER form; dashed line: nuclear PER form.

Figure 1 . 4 :

 14 Figure 1.4: Simulations for the one way switch and toggle switches. Equations are taken from[START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell[END_REF] (see Box 1, Figures 1.e and 1.f in the reference for more details). (A) One way switch: the signal should attain the critical value "Scrit" so that the response passes to the upper state. If the signal decreases, the response stays on the upper state and does not fall back to the lower one. (B) Toggle switch: unlike the one way switch, if the response is on the upper state and the signal decreases, the response falls back to the lower state.

Figure 1 . 5 :

 15 Figure 1.5: Albert Goldbeter cell cycle model: sustained oscillations of cyclin concentration, active cdc2 kinase, and cyclin protease. (B) Limit cycle behavior for the mitotic oscillator obtained by plotting the solution of the system in (cyclin, cdc2) space. Two sets of initial conditions are considered, one inside and the other outside the limit cycle; arrows indicate the direction of the time evolution.

Figure 1 . 6 :

 16 Figure 1.6: Fission yeast cell cycle: circles represent phases of the cell cycle, colored rectangles show the three main transitions and the wiring diagram illustrates the antagonism mechanism between Cdc13/Cdc2 and Rum1, Ste9/APC. Green rectangle: start transition (G1/S). Cdc2/Cdc13 and Ste9/APC, Rum1 mutual inhibition (processes 1,2,3,4); help of the "starter kinase SK" to shift the balance for Cdc2/Cdc13 by deactivating Ste9/APC, Rum1 (processes 5 and 6). Brown rectangle: G2/M transition. Mutual antagonism between Cdc2/Cdc13 and WEE1 (precesses 8 and 9); help of Cdc25 (process 12) by inactivating WEE1 (process 11); Cdc25 is activated in a positive feedback by Cdc2/Cdc13 (process 10). Red rectangle: finish transition (M/G1). Slp1/APC helps Ste9/APC (process 14) by inhibiting the activity of Cdc2/Cdc13 (process 15). Blue and red arrows represent negative feedback loops, which generate the oscillatory activity of Cdc2/Cdc13.

Figure 1 . 7 :

 17 Figure 1.7: Simulations for the cell cycle model proposed by Tyson and Novak: (A) Cell mass (solid line) between birth and division is divided by two when MPF concentration (dashed line) decreases below a threshold value in the end of mitosis. (B) Antagonism between Cdc13 complex (solid line) and its enemies Ste9 (dotted line), and Slp1 (dashed line).

Figure 1 .

 1 Figure 1.8: Simulations for the age structured one-phase model (logarithmic scale). (A) Influence of the transition age on the population growth (logarithmic scale) with and without periodic coefficients. Transition rate K = k 0 1 x>xtrans (A + Bcos(2πt/T )) (B) Comparison of the growth with and without periodic coefficients. Results are coherent with theorem (3), which states that there is no general inequality between λ per and λ s . In particular, one can build examples with λ per < λ s as well as λ per > λ s .

Figure 2 . 2 :Figure 2 . 3 :

 2223 Figure 2.2: Cell cycle dynamics with coupling to the circadian clock. Autonomous period of the cell cycle is equal to 18 h in this example. Increasing the coupling strength tends to regulate the cell cycle to 24 h. Coupling strength: (A) 0, (B) 0.5, (C) 1, (D) 1.5.

Figure 2 . 4 : 6 Figure 2 . 5 :

 24625 Figure 2.4: Effects of coupling the circadian clock to the cell cycle with autonomous period of 20 h. (A) Without coupling (C = 0): MPF activity follows a 20-h autonomous cycle. BMAL1/CLOCK period is equal to 24 h. (B) With coupling (C = 1.2): the cell cycle period is entrained to 24 h. (C) The population in M phase is entrained to 24 h, hence cells have a division cycle of 24 h instead of 20 h. (D) With coupling, the growth rate does not decrease even though the cell cycle period becomes longer.

Figure

  Figure 2.9: Effects of mutating circadian genes. (A) MPF activity under different clock mutations. (B) BMAL1/CLOCK activity under different clock mutations, current model [9]. (C) BMAL1/CLOCK activity under different clock mutations, model of Mirsky et al. [108].

Figure 2 .Figure 2 .

 22 Figure 2.10: (A) Per2 mutants. (B) Wild type cells.

  where x is a cellular dynamical state in an open subset Ω of dimension d and λ a cellular parameter state in an open subset Λ of dimension p. In our case, x = x 1 , • • • , x d and λ represent respectively the set of proteins/mRNAs concentrations and the intrinsic cell cycle period of each cell.

Figure 3 . 1 :

 31 Figure 3.1: Deterministic growth condition.

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: Passive tracer test case. Left figures: particle solution; right figures: upwind-SSP solution. (A-B) Solution at t = 2.5, (C-D) contour plot of the solution at t = 2.5, (E-F) Contour plot of the solution at t = 5.

Figure 3 . 4 :

 34 Figure 3.4: (A) MPF-WEE1 antagonistic activity for λ = 180: division occurs if MPF activity rises above that of WEE1 and then shuts down. (B) Impact of λ on MPF activity. (C) Total number of cells with different source terms with no dependence of MPF activity on the total number of cells (η = 1 × 10 -16 ). Initial number of cells = 4, two cells with λ = 230 and two with λ = 180. (D) Impact of the dependence of MPF activity on the total number of cells (η = 2 × 10 -3 ).

Figure 3 . 5 :

 35 Figure 3.5: Comparison between the average molecular concentrations of all cells and molecular concentrations of a cell chosen randomly in the population (A,B,C) No synchronization between cells: k s = 0; (E,F,G) Synchronized activity: k s = 0.05.

Figure 3 . 6 :

 36 Figure 3.6: (A) Growth rate indicating two peaks of division each day. (B) MPF-WEE1 average activity for all cells. (C) Counting the number of times each cell divides per day. Each square indicates a division. (D) Similar to figure C but simulation shows the division times for 100 cells instead of 10.

Figure A. 1 :

 1 Figure A.1: Parameters estimate for the initial optimization process: one can see how the tumor age is changing with respect to different values of α and K. There exists an interval of values (t, α, K) fitting the biological data.

  Figure A.2: Fitting the absolute growth data using Gompertz model. Solid curve is obtained by finding the appropriate couple values (α, K) yielding an absolute growth (V (t +1)-V (t)) that fit best nakamura data. V (t) is the volume at time t and is given by V c × N (t), where V c is the volume of one cell.

Figure A. 3 :

 3 Figure A.3: Comparison of the absolute growth per year between by Nakamura data and our model. (A) Patient 1, (B) patient 4.

Figure A. 4 :

 4 Figure A.4: First method of fit: evolution of tumor volume in time. From upper left to bottom right: patient 1 till patient 9.

Figure A. 5 :

 5 Figure A.5: Relative growth data fit using Gompertz model.

Figure A. 7 :Figure A. 8 :

 78 Figure A.7: Final approach: evolution of tumor volume in time. From upper left to bottom right: patient 1 till patient 9.

- 8

 8 double i n i t s t a t e v a l u e s f i x e d [ neqs ] = { 0 . 2 , 0 . 4 , 1 . 1 2 , 0 . 8 , 1 . 3 , 0 . 8 , 1 . 1 , 0 . 0 1 5 , 0 . 5 7 5 , 0 . 0 8 } ; double s c a l f a c t 2 f i x [ i n i t c e l l s n u m b e r ] = {

  < b double i n i t s t a t e v a l u e s r a n d [ neqs ][START_REF] Akashi | A positive role for PERIOD in mammalian circadian gene expression[END_REF] = { { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } , { 0 . , 1 . } } ;// Choose t h e l o w e r bound f o r t h e i n t r i n s i n c p e r i o d o f e a ch c e l l : // Chosen randomly b e t w e e n min and max v a l u e s b e l o w :double c e l l c y c l e p e r i o d m i n f = 8 2 . 8 4 ; // 8 2 . 8 4 h o l d s f o r 28 h double c e l l c y c l e p e r i o d m a x f = 1 1 5 . 9 7 ; // 1 1 5 . 9 7 h o l d s f o r 20 h // Parameters f o r t h e mean v a l u e f u n c t i o n double params MVfunction [ 1 0 ] = { 1 , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . } ; #e n d i f / * -

  / #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e <c s t d l i b > #i n c l u d e <cmath> #i n c l u d e < s t d i o . h> #i n c l u d e " haupt . h" void r i g h t h a n d s i d e s a b i n a ( double x , double * * y , double * * dydx , double * params , double * * s c a l f a c t 2 , const i n t neqs , i n t n s t a t e s ) { double y1 , y2 , y3 , y4 , y5 , y6 , y7 , y8 , y9 , y10 , s c a l f a c t o r 2 ; double s c a l f a c t 1 = params [ 0 ] ; double nu1b

  y1 = y [ 0 ] [ j ] ; y2 = y [ 1 ] [ j ] ; y3 = y [ 2 ] [ j ] ; y4 = y [ 3 ] [ j ] ; y5 = y [ 4 ] [ j ] ; y6 = y [ 5 ] [ j ] ; y7 = y [ 6 ] [ j ] ; y8 = y [ 7 ] [ j ] ; y9 = y [ 8 ] [ j ] ; y10 = y [ 9 ] [ j ] ; s c a l f a c t o r 2 = s c a l f a c t 2 [ 0 ] [ j ] ; dydx [ 0 ] [ j ] = nu1b * ( y7 + c ) / ( k1b * (1+pow ( y3 / k 1 i , p ) ) + y7 + c ) -k1d * y1 + MVfunction out [ 0 ]

  / // Param f o r t h e c i r c a d i a n c l o c k ( s o u r c e : Becker-weimann e t a l . ) and // t h e c e l l c y c l e ( El C h e i k h e t a l .2 0 1 4 ) #i n c l u d e " haupt . h" #i n c l u d e <cmath> #i n c l u d e <math . h> #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e < s t d i o . h> void O D E s y s t e m p a r a m s a s s i g n a t i o n ( double * params ) { /

  2 4 / params [ 0 ] ; // k 2 t params [ 1 1 ] = . 0 2 / params [ 0 ] ; // k 3 t params [ 1 2 ] = . 1 2 / params [ 0 ] ; // k3d params [ 1 3 ] = . 6 / params [ 0 ] ; // nu4b params [ 1 4 ] = . 1 6 ; // k 4 b params [ 1 5 ] = . ; // r params [ 1 6 ] = . 7 5 / params [ 0 ] ; // k4d params [ 1 7 ] = . 2 4 / params [ 0 ] ; // k 5 b params [ 1 8 ] = . 0 6 / params [ 0 ] ; // k5d params [ 1 9 ] = . 4 5 / params [ 0 ] ; // k 5 t params [ 2 0 ] = . 0 6 / params [ 0 ] ; // k 6 t params [ 2 1 ] = . 1 2 / params [ 0 ] ; // k6d params [ 2 2 ] = . 0 9 / params [ 0 ] ; // k6a params [ 2 3 ] = . 0 0 3 / params [ 0 ] ; // k7a params [ 2 4 ] = . 0 9 / params [ 0 ] ; // k7d // Param f o r t h e mpf-wee1 d y n a m i c s params [ 2 5 ] = ; // nmw params [ 2 6 ] = . ; // k0mpf : 10 = 4 + 6 * e x p (-lambda * N) params [ 2 7 ] = . 0 5 ; // k1mpf params [ 2 8 ] = . 5 ; // k1wee1 params [ 2 9 ] = 5 . 0 ; // dwee1 params [ 3 0 ] = 1 . 0 ; // k a c t w params [ 3 1 ] = 2 0 0 . 0 ; // k i n a c t w params [ 3 2 ] = 1 . ; // dw1 params [ 3 3 ] = 1 . ; // dw2 params [ 3 4 ] = 1 . 2 ; // kbm C o u p l i n g f o r c e params [ 3 5 ] = 0 . 5 ; // cw params [ 3 6 ] = 0 . 0 1 ; // k a c t r a t e c o n s t a n t f o r i n h i b i t o r params [ 3 7 ] = 5 0 ; // s s t r e n g t h o f t h e i n h i b i t o r ; params [ 3 8 ] = 1 . e -16; // lambda k0mpf = 4 + 6 * e x p (-lambda * n s t a t e s ) } / * -

  -RK4 step . cpp

  / // Runge k u t t a f u n c t i o n o v e r a t i m e s t e p t a u#i n c l u d e " haupt . h" #i n c l u d e <cmath> #i n c l u d e <math . h> #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e < s t d i o . h> u s i n g s t d : : cou t ; u s i n g s t d : : e n d l ; u s i n g s t d : : setw ; void RK4 step ( const double tau , double x , double * * y ,double * * dydx , double * params , double * * s c a l f a c t 2 , const i n t neqs , i n t n s t a t e s ) { double * * k1 , * * k2 , * * k3 , * * k4 ; double * * y aux ;y aux = md2 new ( neqs , n s t a t e s ) ; k1 = md2 new ( neqs , n s t a t e s ) ; k2 = md2 new ( neqs , n s t a t e s ) ; k3 = md2 new ( neqs , n s t a t e s ) ; k4 = md2 new ( neqs , n s t a t e s ) ;f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < neqs ; i++ ) { y aux [ i ] [ j ]= y [ i ] [ j ] ; } } r i g h t h a n d s i d e s a b i n a ( x , &y aux [ 0 ] , &dydx [ 0 ] , &params [ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , n s t a t e s ) ; f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < neqs ; i++ ) { k1 [ i ] [ j ] = dydx [ i ] [ j ] ; y aux [ i ] [ j ] = y [ i ] [ j ] + 0.5 * tau * k1 [ i ] [ j ] ; } } r i g h t ha n d s i d e s a b i n a ( x + 0 . 5 * tau , &y aux [ 0 ] , &dydx [ 0 ] , &params [ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , n s t a t e s ) ;f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < neqs ; i++ ) { k2 [ i ] [ j ] = dydx [ i ] [ j ] ; y aux [ i ] [ j ] = y [ i ] [ j ] + 0.5 * tau * k2 [ i ] [ j ] ; } }r i g h t h a n d s i d e s a b i n a ( x + 0 . 5 * tau , &y aux [ 0 ] , &dydx [ 0 ] , &params [ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , n s t a t e s ) ;f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < neqs ; i++ ) { k3 [ i ] [ j ] = dydx [ i ] [ j ] ; y aux [ i ] [ j ] = y [ i ] [ j ] + tau * k3 [ i ] [ j ] ; } } r i g h t ha n d s i d e s a b i n a ( x + tau , &y aux [ 0 ] , &dydx [ 0 ] , &params [ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , n s t a t e s ) ;f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < neqs ; i++ ) { y[ i ] [ j ] += ( tau * k1 [ i ] [ j ] + 2 . * tau * k2 [ i ] [ j ] + 2 . * tau * k3 [ i ] [ j ] + tau * dydx [ i ] [ j ] ) / 6 . ; } } m d 2 fr e e ( y aux , neqs , n s t a t e s ) ; m d 2 f r e e ( k1 , neqs , n s t a t e s ) ; m d 2 f r e e ( k2 , neqs , n s t a t e s ) ; m d 2 f r e e ( k3 , neqs , n s t a t e s ) ; m d 2 f r e e ( k4 , neqs , n s t a t e s ) ; } / * -

  main prgram . cpp

  / #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e < s t d l i b . h> #i n c l u d e <f s t r e a m > #i n c l u d e <cmath> #i n c l u d e <c s t r i n g > #i n c l u d e <time . h> #i n c l u d e < s t d i o . h> #i n c l u d e <deque> #i n c l u d e <random> #i n c l u d e " haupt . h" #i n c l u d e " p a r a m e t e r s . h" u s i n g namespace s t d ; i n t main ( ) { i n t n s t a t e s , n n s t a t e s , a d d e d s i z e , Npop , N p o p t o t a l ; const double tau = t i m e s t e p ; n s t a t e s = i n i t c e l l s n u m b e r ;

s c a l f a c t 2 =

 2 t h m double * params ; // t r a c k s t h e number o f d i v i d e d c e l l s deque< i n t > d i v i d e d c e l l s n u m b e r ; deque< i n t > d i v i d e d c e l l s n u m b e r o u t p u t ; // t r a c k s t h e t i m e each c e l l d i v i d e d deque< double > d i v i d e d c e l l s t i m e ; deque< double > d i v i d e d c e l l s t i m e o u t p u t ; deque <double> KG f1 ; // d-v e c t o r t h a t c a c l c u l a t e s m o l e c u l a r a v e r a g e s a t each t i m e s t e p deque <double>MVfunction output ; md2 new ( 1 , n s t a t e s ) ; d i v c o u n t e r 1 = mi2 new ( 1 , n s t a t e s ) ; d i v c o u n t e r 2 = mi2 new ( 1 , n s t a t e s ) ; params = md new ( nparams ) ; K = md2 new ( neqs , n s t a t e s ) ; // Each column o f K r e p r e s e n t s t h e m o l e c u l a r s t a t e o f a c e l l dydx = md2 new ( neqs , n s t a t e s ) ; // R i g h t hand s i d e o f ODE s y s t e m // I n i t i a t i o n i n t i n d 1 = 0 ; i n t i n d 2 = 0 ; i n t i n d d i v ; a d d e d s i z e = 0 ; Npop = 0 ; N p o p t o t a l = n s t a t e s ;

  c l o s e ( ) ; o f s t r e a m w r i t e t i m e ( " time . t e x t " ) ; o f s t r e a m w r i t e o n e c e l l ( " o n e c e l l . t e x t " ) ; o f s t r e a m w r i t e c e l l s n u m b e r ( " c e l l s n u m b e r . t e x t " ) ; o f s t r e a m w r i t e m o l e c u l a r a v e r a g e s ( " m o l e c u l a r a v e r a g e s . t e x t " ) ; o f s t r e a m o n e c o m p d i f f c e l l s ( " o n e c o m p d i f f c e l l s . t e x t " ) ; i n t c o u n t e r = 0 ; double Time = 0 ; while ( Time < FTimetau && n n s t a t e s < 2 0 0 0 0 ) { w r i t e t i m e << Time << e n d l ; Time+=t a u ; c o u n t e r +=1;

  Runga-Kutta s o l v e r RK4 step ( tau , Time , &K[ 0 ] , &dydx [ 0 ] , &params [ 0 ] , &s c a l f a c t 2 [ 0 ] , neqs , n n s t a t e s ) ;// Examines i f d i v i s i o n c o n d i t i o n s a r e s a t i s f i e d d i v i s i o n t y p e ( d i v i s i o n t y p e , Time , tau, K G f u n c t i o n s t i f f , K G d i v i s i o n c o n s t a n t , index MPF , index WEE1 , &K[ 0 ] , &d i v c o u n t er 1 [ 0 ] , &d i v c o u n t e r 2 [ 0 ] , n n s t a t e s , Npop , d i v i d e d c e l l s n u m b e r , d i v i d e d c e l l s n u m b e r o u t p u t , d i v i d e d c e l l s t i m e , d i v i d e d c e l l s t i m e o u t p u t , KG f1 ) ; // c o u t << KG f1 [ c o u n t e r -1] << e n d l ; o n e c o m p d i f f c e l l s << K [ 6 ] [ 0 ] << " " << K [ 6 ] [ 1 ] << " " << K [ 6 ] [ 2 ] d i v i s i o n t y p e . cpp

  / #i n c l u d e " haupt . h" #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e < s t d l i b . h> #i n c l u d e <f s t r e a m > #i n c l u d e <cmath> #i n c l u d e <c s t r i n g > #i n c l u d e <time . h> #i n c l u d e < s t d i o . h> #i n c l u d e <random> #i n c l u d e <deque> u s i n g namespace s t d ; void d i v i s i o n t y p e ( char * d i v i s i o n t y p e , double Time , const double tau , const double s t i f f n e s s , const double t r , i n t index MPF , i n t index WEE1 , double * * K, i n t * * d i v c o u n t e r 1 , i n t * * d i v c o u n t e r 2 , i n t n n s t a t e s , i n t &Npop , deque <int> &d i v i d e d c e l l s n u m b e r , deque <int> &d i v i d e d c e l l s n u m b e r o u t p u t , deque <double> &d i v i d e d c e l l s t i m e , deque <double> &d i v i d e d c e l l s t i m e o u t p u t , deque <double> &KG f1 ) { double r ; // random number fo r t h e k o s h l a n d d i v i s i o n t y p e / * i n i t i a l i z e random s e e d : * / // s r a n d ( t i m e (NULL ) ); i f ( strcmp ( d i v i s i o n t y p e , " d i r a c " ) == 0) { KG f1 . push back ( f u n c t i o n G ( 0 . 0 6 , K[ index MPF ] [ 2 ] , s t i f f n e s s , s t i f f n e s s ) ) ; f o r ( i n t j = 0; j < n n s t a t e s ; j ++) { i f ( K[ index MPF ] [ j ] > K[ index WEE1 ] [ j ] ) { d i v c o u n t e r 1 [ 0 ] [ j ] += 1 ; d i v c o u n t e r 2 [ 0 ] [ j ] = 0 ; // i n d d i v = 0 ;// c o u t << d i v c o u n t e r 1 [ 0 ] [ j ] << e n d l ;} i f ( (K[ index MPF ] [ j ] < K[ index WEE1 ] [ j ] ) && (K [ index MPF ] [ j ] < 0 . 0 6 ) && ( d i v c o u n t e r 1 [ 0 ] [ j ] >= 1) ){ d i v c o u n t e r 2 [ 0 ] [ j ] += 1 ; } i f ( d i v c o u n t e r 2 [ 0 ] [ j ] == 1) { Npop +=1;// c o u t << Npop << s e t w ( 1 2 ) << Time << e n d l ; d i v i d e d c e l l s n u m b e r . push back ( j ) ; d i v i d e d c e l l s n u m b e r o u t p u t . push back ( j ) ; d i v i d e d c e l l s t i m e . push back ( Time ) ; d i v i d e d c e l l s t i m e o u t p u t . push back ( Time ) ; // c o u t << Npop << s e t w ( 1 2 ) << d i v i d e d c e l l s n u m b e r . s i z e ( ) << e n d l ; d i v c o u n t e r 1 [ 0 ] [ j ] = 0 ; d i v c o u n t e r 2 [ 0 ] [ j ] = 0 ; } // c o u t << d i v c o u n t e r 2 [ 0 ] [ j ] << s e t w ( 1 2 ) << i n d d i v << s e t w ( 1 2 ) << Npop << e n d l ; } } e l s e i f ( strcmp ( d i v i s i o n t y p e , " koshland " ) == 0) { } e l s e { c o u t << " type not d e f i n e d " << e n d l ; } } / * -

  / // F u n c t i o n t h a t u p d a t e t h e s t a t e o f t h e s y s t e m and a l l o c a t e new memory // due t o p o s s i b l e i n c r e a s e i n t h e c e l l number#i n c l u d e " haupt . h" #i n c l u d e <i o s t r e a m > #i n c l u d e <iomanip> #i n c l u d e < s t d l i b . h> #i n c l u d e <f s t r e a m > #i n c l u d e <cmath> #i n c l u d e <c s t r i n g > #i n c l u d e <time . h> #i n c l u d e < s t d i o . h> #i n c l u d e <deque> u s i n g namespace s t d ;void s t a t e s u p d a t e ( double Time , char * i n i t i a l v a l u e s , double i n i t s t a t e v a l u e s f [ ] , const i n t neqs , i n t &a d d e d s i z e , i n t &n s t a t e s , i n t &n n s t a t e s ,i n t &i n d d i v , i nt &N p o p t o t a l , i n t &Npop , double * * K, double * * dydx , double * * s c a l f a c t 2 , i n t * * d i v c o u n t e r 1 , i n t * * d i v c o u n t e r 2 , deque <int> &d i v i d e d c e l l s n u m b e r , deque <double> &d i v i d e d c e l l s t i m e , double c e l l c y c l e p e r i o d m i n f , double c e l l c y c l e p e r i o d m a x f , o f s t r e a m &w r i t e c e l l s n u m b e r ) {

  a d d e d s i z e = Npop ; i f ( ( a d d e d s i z e > 0 ) ) { K = m d 2 r e a l l o c (&K[ 0 ] , neqs , n n s t a t e s , 0 , a d d e d s i z e ) ; dydx = m d 2 r e a l l o c (&dydx [ 0 ] , neqs , n n s t a t e s , 0 , a d d e d s i z e ) ; d i v c o u n t e r 1 = m i 2 r e a l l o c (& d i v c o u n t e r 1 [ 0 ] , 1 , n n s t a t e s , 0 , a d d e d s i z e ) ; d i v c o u n t e r 2 = m i 2 r e a l l o c (& d i v c o u n t e r 2 [ 0 ] , 1 , n n s t a t e s , 0 , a d d e d s i z e ) ; s c a l f a c t 2 = m d 2 r e a l l o c (& s c a l f a c t 2 [ 0 ] , 1 , n n s t a t e s , 0 , a d d e d s i z e ) ; // Adding t h e new c e l l s : n n s t a t e s += a d d e d s i z e ; N p o p t o t a l +=Npop ;

  [ i ] [ j ] = K[ i ] [ d i v i d e d c e l l s n u m b e r [ jn s t a t e s ] ] ; // K[ i ] [ j ] = i n i t s t a t e v a l u e s f [ i ] + r a n d o m f l o a t ( 0 . , 0 . 2 ) ; // c o u t << << e n d l ;} s c a l f a c t 2 [ 0 ] [ j ] = r a n d o m f l o a t ( c e l l c y c l e p e r i o d m i n f , cell c y c l e p e r i o d m a x f ) ; d i v c o u n t e r 1 [ 0 ] [ j ] = d i v c o u n t e r 1 [ 0 ] [ d i v i d e d c e l l s n u m b e r [ jn s t a t e s ] ] ; d i v c o u n t e r 2 [ 0 ] [ j ] = d i v c o u n t e r 2 [ 0 ] [ d i v i d e d c e l l s n u m b e r [ jn s t a t e s ] ] ; } i f ( ptr == 0) { c o u t << " no s p a c e s h o u l d be a l l o c a t e d b e c a u s e n i s NULL \n" ; } return ( ptr ) ; } void m d f r e e ( double * a , i n t n ) { f r e e ( a ) ; } // A l l o c a t i o n and f r e e i n g 2D double * * md2 new ( i n t m, i n t n ) { double * * a = new double * [m ] ;f o r ( i n t i = 0; i < m; i ++) a [ i ] = new double [ n ] ; return a ; } i n t * * mi2 new ( i n t m, i n t n ) { i n t * * a = new i n t * [m ] ; f o r ( i n t i = 0; i < m; i ++) a [ i ] = new i n t [ n ] ; return a ; } void m d 2 f r e e ( double * * a , i n t m, i n t n ){ f o r ( i n t i = 0; i < m; i ++) d e l e t e [ ] a [ i ] ; d e l e t e [ ] a ; } void m i 2 f r e e ( i n t * * a , i n t m, i n t n ){ f o r ( i n t i = 0; i < m; i ++) d e l e te [ ] a [ i ] ; d e l e t e [ ] a ; } // c a l l o c double * * m d 2 c a l l o c ( i n t m, i n t n ) { double * * xyz ; xyz = ( double * * ) c a l l o c (m, s i z e o f ( double * ) ) ; f o r ( i n t i = 0; i < m; i ++) xyz [ i ] = ( double * ) c a l l o c (n , s i z e o f ( double ) ) ; return xyz ; return xyz ; } // A l l o c a t i o n and f r e e i n g 3D // a l l o c a t e a 3D a r r a y double * * * md3 new ( i n t x , i n t y , i n t z ) { double * * * t h e a r r a y = new double * * [ x ] ; f o r ( i n t i ( 0 ) ; i < x ; ++i ) { t h e array [ i ] = new double * [ y ] ; f o r ( i n t j ( 0 ) ; j < y ; ++j ) { t h e a r r a y [ i ] [ j ] = new double [ z ] ; f o r ( i n t k ( 0 ) ; k < z ; ++k ) { t h e a r r a y [ i ] [ j ] [ k]= 0 . ; } 144 } } return t h e a r r a y ; } void m d 3 f r e e ( double * * * t h e a r r a y , i n t x , i n t y , i n t z ) { f o r ( i n t i = 0; i < x ; ++i ) { f o r ( i n t j = 0; j < y ; ++j ) { d e l e t e [ ] t h e a r r a y [ i ] [ j ] ; } d e l e t e [ ] t h e a r r a y [ i ] ; } d e l e t e [ ] t h e a r r a y ; }// R e a l l o c a t i o n o f a m a t r i x o f s i z e [m ] [ n ] // New s i z e i s [m + mm] [ n + nn] double * * m d 2 r e a l l o c ( double * * xyz , i n t m, i n t n , i n t mm, i n t nn ) { xyz = ( double * * ) r e a l l o c ( xyz , (m + mm) * s i z e o f ( double * ) ) ; // The new column ' s p o i n t e r must be i n i t i a l i s e d t o NULL f o r ( i n t i = m; i < (m+mm) ; i ++) xyz [ i ] = NULL ; // R e a l l o c a t e rows f o r ( i n t i = 0; i < (m+mm) ; i ++) // xyz [ i ] = ( i n t * ) r e a l l o c ( xyz [ i ] , (n + nn ) * s i z e o f ( i n t ) ) ; xyz [ i ] = ( double * ) r e a l l o c ( xyz [ i ] , (n + nn) * s i z e o f ( double ) ) ; return xyz ; } // I n t e g e r m a t r i x r e a l l o c a t i o n i n t * * m i 2 r e a l l o c ( i n t * * xyz , i n t m, i n t n , i n t mm, i n t nn ) { xyz = ( i n t * * ) r e a l l o c ( xyz , (m + mm) * s i z e o f ( i n t * ) ) ; // The new column ' s p o i n t e r must be i n i t i a l i s e d t o NULL f o r ( i n t i = m; i < (m+mm) ; i ++) xyz [ i ] = NULL ; // R e a l l o c a t e rows f o r ( i n t i = 0; i < (m+mm) ; i ++) // xyz [ i ] = ( i n t * ) r e a l l o c ( xyz [ i ] , (n + nn ) * s i z e o f ( i n t ) ) ; xyz [ i ] = ( i n t * ) r e a l l o c ( xyz [ i ] , (n + nn) * s i z e o f ( i n t ) ) ; return xyz ; } Multiscale modeling for the regulation of the cell cycle by the circadian clock : Applications to chronotherapy.

  Image en couverture : LU : Arnold tongues. RU : Devil's staircase. LD : Growth rate. RD : divisions count per day

  and S. Bernard. Modeling Biological Rhythms in Cell Populations. Mathematical Modelling of Natural Phenomena / Volume 7 / Issue 06 / 2012, pp 107 -125. • Luc Paquet, Raouf El Cheikh, Dominique Lochegnies and Norbert Siedow. Radiative Heating of a Glass Plate. MathematicS in Actions. Vol.5 no.1 (2012), p. 1-30.

  . . . . . 2.2 Cell cycle dynamics with and without coupling to the circadian clock . . . . 2.3 Arnold tongues and devil's staircase . . . . . . . . . . . . . . . . . . . . . . . 2.4 Effects of coupling the circadian clock to the cell cycle on MPF-WEE1 dynamics, on the fraction of cells in mitosis and on the Growth rate. Autonomous cell cycle period is equal to 20 h . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Effects of coupling the circadian clock to the cell cycle on the growth rate . . 2.6 Effects of mutating circadian genes on the fraction of cells entering mitosis . 2.7 Effects of mutating circadian genes on the growth rate . . . . . . . . . . . . 2.8 Effects of mutating circadian genes on the growth rate. Simulations with other circadian models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Effects of mutating circadian genes on MPF and BMAL1/CLOCK activity . 2.10 Comparison of transition rates between Per2 mutants and wild type cells . . 2.11 Comparison of transition rates between Cry2 mutants and wild type cells . . 3.1 Deterministic growth condition. . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Passive tracer numerical test case for the particle method . . . . . . . . . . . 3.3 Solution of the molecular structured equation without a source term . . . . . 3.4 Total number of cells with different source terms . . . . . . . . . . . . . . . . 3.5 Synchronization: influence of the dependence of the rate of Per2/Cry production on the its average concentration . . . . . . . . . . . . . . . . . . . . . . 3.6 Influence of the parameter λ on the growth rate . . . . . . . . . . . . . . . . A.1 Parameters estimate for the initial optimization process . . . . . . . . . . . . A.2 Absolute Growth fit using Gompertz model . . . . . . . . . . . . . . . . . . A.3 Comparison of the absolute growth per year between Nakamura data and our model for two patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 First method of fit: evolution of tumor volume in time . . . . . . . . . . . . A.5 Relative growth fit using Gompertz model . . . . . . . . . . . . . . . . . . . A.6 Second method of fit: evolution of tumor volume in time . . . . . . . . . . . A.7 Combined approach: evolution of tumor volume in time . . . . . . . . . . . . A.8 Prediction of the tumor growth in time . . . . . . . . . . . . . . . . . . . . .

  The model of Leloup et al. includes only Bmal1, Per, Cry and Rev-erbα. A good test for the validity of the model is its ability to reproduce correct phases of circadian variables. This model captures precise phase relationships among molecular components of the circadian clock work which is not always the case for other models. Namely: Rev-erbα mRNA leads Per1, and Per2 mRNAs by 4 h; Per1, and Per2 mRNAs lead Cry1, Cry2, and Rorc mRNAs by 4 h; Cry1, Cry2, and Rorc mRNAs lead Clk and Bmal1 mRNAs by 4 h; Clk and Bmal1 mRNAs lead Rev-erbα mRNA by 8 h. The reason is that the authors made exclusive use of cell-level data for both parametric fitting and validation which make it more accurate at the cell-level where other models fail.

	al.
	and Forger et al. this model has several advantages. It includes eight genes: Per1, Per2,
	Cry1, Cry2, Clk, Bmal1, Rev-erbα and Rorc, whereas the model of Forger et al. includes
	Per1, Per2, Cry1, Cry2, Rev-erbα and BMAL1-CLK which are implicitly present at high
	constant levels.

  This result needs explanation. Suppose we have two sets of coefficients d 1 , d 2 , B 1 , B 2 with the same period T , and the corresponding growth rates λ 1 per , λ 2 per . We define an intermediate model with coefficients d θ , B θ , with θ ∈ (0, 1), by

.4.21) λ per is geometrically convex with respect to B and convex with respect to d.

  The dynamical variables of the circadian clock are: y 1 Per2 or Cry mRNA and proteins; y 2 PER2/CRY complex (cytoplasm); y 3 PER2/CRY complex (nucleus); y 4 Bmal1 mRNA; y 5 BMAL1 cytoplasmic protein; y 6 BMAL1 nuclear protein; y 7 Active BMAL1; The dynamical variables of the cell cycle are: z 8 Active MPF; z 9 Active WEE1; z 10 Active MPF inhibitor.

						.2.6)
	dy 7 dt	= k 6a y 6 -k 7a y 7 -k 7d y 7 ,	(2.2.7)
	dz 8 dt	=	k 0mpf k n 1mpf 1mpf + z n k n 8 + sz n 10	(1 -z 8 ) -d wee1 z 9 z 8 ,	(2.2.8)
	dz 9 dt	=	k actw k actw + d w1	(c w + Cy 7 )+
		k actw k actw + d w1	-1	k inactw z n 8 z 9 k n 8 1wee1 + z n	-d w2 z 9 ,	(2.2.9)
	dz 10 dt	= k act (z 8 -z 10 ).	(2.2.10)

Table 2 .

 2 1: Effects of mutations on the period of the circadian clock: comparison between experimental data and simulations.

	Mutation	Circadian clock period
		Experimental	Simulation
	Per2	Arrhythmic [108]	Arrhythmic
	Bmal1	Arrhythmic [108]	Arrhythmic
	Cry2	Rhythmic, long period [108] Rhythmic, T = 24.2 h
	Per2 /Cry2 Rhythmic [119]	Rhythmic, T = 22.7 h

Table 2 .

 2 Figure 2.6: Effects of mutating circadian genes on the fraction of cells entering mitosis. (A) Autonomous period of the cell cycle is equal to 28 h. (B) Autonomous period of the cell cycle is equal to 20 h. 2: Effects of mutations on the period of the M phase (column period) and the growth rate (column g.r.).

		1.6						1.6
					WT									WT
	Fraction of cells in M phase	0.2 0.4 0.6 0.8 1.0 1.2 1.4			mCry2 mCry1/mCry2 mBmal1 mPer2/mCry2 mPer2	Fraction of cells in M phase	0.2 0.4 0.6 0.8 1.0 1.2 1.4						mCry2 mCry1/mCry2 mBmal1 mPer2/mCry2 mPer2
		0.0						0.0
		1.5	2	2.5	3	3.5	4	0	0.5	1	1.5	2	2.5	3	3.5	4
			Days				(B)					Days
			Autonomous period 28 h Autonomous period 20 h
		mutation	period (h)	g.r. (d -1 ) period (h)		g.r. (d -1 )
		Wild-type			28.0		0.3389		20.0		0.3282
		Per2			56.4		0.3364		38.4		0.3954
		Bmal1			60.0		0.2938		42.7		0.3558
		Cry2			24.2		0.1895		24.2		0.0986
		Per2 /Cry2			22.6		0.3348		22.7		0.3289

  al.[START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF] and Leloup et al.[START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF]. Results given by the model proposed by Mirsky et al.[START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF] were consistent with current model. Namely, Per2 mutation increases the growth rate, Cry1 decreases it, Bmal1 mutation increases it for autonomous periods less than 22 h and decreases it elsewhere. Simulations on Per and Cry mutations done with the model proposed byLeloup et al.

		0.8
							WT
		0.7					mPer2 case1
							mPer2 case2
		0.6					mBmal1
							mCry2
	Growth rate	0.2 0.3 0.4 0.5					mPer2/mCry2
		0.1
		0.0
		-0.1
		0	10	20	30	40	50	60
			Autonomous period of the cell cycle
					(hours)
					68

Figure

2

.7: Effects of mutating circadian genes on the growth rate. Mutation of Per2, assuming its role in the negative feedback loop (dot-dashed line), decreases the growth rate for periods ranging from 27 to 31 h, for periods larger than 40 h and increases it elsewhere. Mutation of Per2 assuming its positive regulation of Bmal1 (grey dot-dashed line) increases the growth rate almost everywhere. Bmal1 knockout (dotted line) increases the growth rate for periods shorter than 21 h and decreases it elsewhere. Cry2 mutation (grey dashed line) decreases the growth rate everywhere. Per2 /Cry2 double mutation (dashed line) maintains a normal proliferation. et

  .11).

		0.8
		0.5								WT				0.7 0.7
		0.4 0.7 0.6								WT mPer2 mCry1 mCry2 mBmal1 mCry2 mCry1/mCry2 mBmal1 mPer2 mPer2/mCry2		0.6						WT mBmal1 mCry mPer
	Growth rate MPF	0.3 0.3 0.5 0.4											Growth rate	0.5
		0.2 0.2												0.4
		0.1
		0.1 0.0												0.3
		0 0	10 0.5	1	20	1.5	30	2	40 2.5	50 3	3.5	60	4	0	10	20	30	40	50	60
	(A) (A)		Autonomous period of the cell cycle (hours) Days	(B) (B)	Autonomous period of the cell cycle (hours)
			Mutation					Circadian clock rhythmicity
									Current				Mirsky et al.
			Per2				Arrhythmic		Arrhythmic
			Bmal1				Arrhythmic		Arrhythmic
			Cry2				Rhythmic, T = 24.2 h Rhythmic, T = 32.1 h
			Per2 /Cry2 Rhythmic, T = 22.7 h Arrythmic
			Cry1				Arrhythmic		Arryhthmic
	Table 2.3: Effects of mutations on the period of the circadian clock: comparison between
	the current model [9] and the model of Mirsky et al. [108].

Figure

2

.8: Effects of mutating circadian genes on the growth rate. Simulations were performed with other models chosen from literature for the circadian clock (Mirsky et al.

[START_REF] Mirsky | A model of the cell-autonomous mammalian circadian clock[END_REF]

(A), and Leloup and Goldbeter

[START_REF] Leloup | Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF]

(B)). Both models gave a result for Bmal1 mutants similar to our result, Bmal1 mutation increases the growth rate for autonomous periods less than 22 hours and decreases it elsewhere. For other types of mutations, the first model is more consistent with our model. It gave similar result for Per2 mutants, Per2 mutation increases the growth rate (dot-dashed line (A)). Results for Cry1 mutation were also in agreement with our results. Cry1 mutation (grey solid line (A)) seems to decrease the growth rate for a large interval of autonomous periods. Results for Cry2 mutation are not in agreement with ours, since Cry2 mutation (grey dashed line (A)) seems to increase the growth rate and do not decrease it as it was predicted by our model. Simulations on Per and Cry mutations done with the model proposed by Leloup et al. did not show much difference for the growth rate compared to wild type cells (dot-dashed and grey solid lines (B)).

  al. were more consistent with our results. The model proposed by Leloup et al. did not show a difference between mutated and wild type cells. This may be explained by the fact that the model proposed by Mirsky et al. was designed to study the effects of mutating circadian genes, whereas the model proposed by Leloup et al. was designed to generate sustained oscillations, which makes it more robust about parameter variation.

  Bmal1 mutation results in an arrhythmic clock with zero concentration of BMAL1/CLOCK and Per2 mutation results in an arrhythmic clock with a low BMAL1/CLOCK concentration (Figures 2.9B,C and Tables 2.3, 2.4). Consequently, Bmal1 and Per2 mutations have similar effects on the growth rates (Figures 2.7 and 2.8A, dotted and dash-dotted lines). Cry2 mutation results in a rhythmic clock with a longer period for both models. However, the period obtained with the model of Mirsky et al. is longer than with the current model (T = 32.1 h vs T = 24.2 h, Table 2.3). The effect on BMAL1/CLOCK concentration is different: the Mirsky et al. model showed a lower concentration compared to wild type, while the current model showed a higher concentration (Figures 2.9B,C and Table 2.4). This results in different growth rates for the two models (Figures 2.7 and 2.8A, grey dashed lines). Per2 /Cry2 double mutation results in a rhythmic clock with the current model (T = 22.7 h), with a concentration similar to wild type, while it results in an arrhythmic clock using the model of Mirsky et al. (Figures 2.9B,C and Tables 2.3, 2.4). Finally, both models gave an arrhythmic clock for Cry1 mutation, with higher BMAL1/CLOCK concentration compared to wild type (Figures 2.9B,C, grey solid line).

Table A .

 A 1: Biological data: The first column gives the birth date of each patient and the second one gives the date of operation. The average age of cells in the third column was calculated using carbon14 C dating. The growth rate given by the Ki-67 is equal to Ki-67 index cycle time ; if the Ki-67 index is equal to 0.3 and the length of the cycle is 2 days for example, then the growth rate per year is equal to 0.3 2×100 × 365. The last column gives the volume of the tumor when it was extracted during the operation. boundary condition n(0, t) accounts for the addition of new cells in time following a Gompertz type of growth. N (t) is the total number of cells and is obtained by summing up all cells of all ages, N (t) =

	Patient	Birthday		Date of operation	Average age of cells	Ki-67 index	Tumor size(cm3)
	1		23.Aug.-1932	07.Apr.-2004	2.5	1.00	50.5
	2		26.May-1933	09.Jan.-2001	3.3	0.30	1.7
	3		31.Dec.-1939	09.Nov.-2001	4	1.1	44.5
	4		26.Dec.-1948	04.Dec.-1996	3.7	0.7	7.5
	5		07.July.-1953	08.Jan.-2003	1.7	0.5	6.25
	6		02.Sep.-1963	15.Jan.-2002	2.5	0.4	2.25
	7		08.Oct.-1966	20.June.-2011	6	0.7	38.5
	8		03.Jan.-1969	15.Feb.-2011	5	0.5	21
	9		09.June.-1967	21.Mar.-2011	0.5	7.10	168
		t
			n(x, t)dx. The average age of cells can be calculated using the following
	formula a =	0 0 xn(x, t)dx t N (t)	.

  .2, Figures A.4, A.8).

  Figure A.6: Second method of fit: evolution of tumor volume in time. From upper left to bottom right: patient 1 till patient 9.even though it has a small size (see for example patients 2 and 6 in Figure A.7). We finally give a prediction for the evolution of tumor over 25 years, imagining that no surgery was done(Figure A.8). By examining whether the growth curve is in its steady state or exponential state, practitioners could use our model to complement medical imagery about the decision of surgery time.Table A.2: Simulation results using the first method of fit (estimate of the absolute growth).TableA.3: Simulation results using the second method of fit (estimate of the relative growth).

	Patient Estimated proliferation coefficient α Estimated carrying capacity K Tumor age
												First fit
					0.17				346.4707		16.8419
					0.0750				58.4924		27.2273
					0.1410				177.8720		22.5013
					0.0940				77.3119		26.3729
					0.1480				186.4637		14.2414
					0.1140				99.9973		17.2317
					0.1160				90.9549		29.8443
					0.1040				90.5833		29.1302
					0.4830				2.6777e+06		2.7450
	Patient Estimated proliferation coefficient α Estimated carrying capacity K Tumor age
												Second fit
					1.4043				70.3386		5.1446
					0.0786				78.8727		26.6122
					0.6379				63.3347		9.6569
					0.1961 0.9663 0.4888 0.2832 0.2194				19.8495 9.3670 4.5068 47.9088 Tumor size (cm3) 38.6553		17.6683 5.1240 7.9667 17.9972 19.1138
					0.0145				6.7521		0.2001
	P	Experimental data	Nakamura growth rates		Fitted data		Parameters estimate
		< a >	V	Ki-67	γr	γa	Ṽ	< ã >	γr	γa	α	K	t
	1	2.5000	50.5000	1,00	2.5554	0.1682	50.6392	2.4773	1.7582	0.3013	1.7775	69.7568	4.5906
	2	3.3000	1.7000	0,30	23.8801	0.2415	1.6935	3.2964	20.9658	0.3698	0.1250	12.5342	21.0108
	3	4.0000	44.5000	1,10	3.3508	0.4100	44.6559	3.9713	2.4459	0.6760	0.7872	61.9903	8.6434
	4	3.7000	7.5000	0,70	14.5470	1.1516	7.4916	3.7070	14.5561	1.0613	0.1943	19.9917	17.7718
	5	1.7000	6.2500	0,50	15.6935	1.0738	6.1374	1.7485	17.9387	0.8890	0.7610	10.0796	5.9235
	6	2.0000	2.2500	0,40	22.1176	1.9479	2.3334	2.2154	33.1126	0.8155	0.1649	21.3576	15.1021
	7	6.0000	38.5000	0,70	4.2615	1.4842	38.4998	6.0001	4.3330	1.4833	0.2795	48.1271	18.1287
	8	5.0000	21.0000	0,50	8.0728	1.6546	20.9990	5.0005	8.3744	1.6480	0.2097	39.8848	19.6095
	9	0.5000	168.0000	7,10	-5.0026	1.5706	167.8469	0.7581	4.2916	1.4473	5.5687	215.6970	1.4394
						Table A.4: Final results
							121

Table B .

 B 2: Parameter values for Tyson and Novak fission yeast cell cycle model. All parameters have units min -1 , except the Ji s and k diss which are dimensionless constants. Table B.3: Parameters description for the coupled population-molecular model (Section 2.2).

		const i n t neqs = 1 0 ;
		// Number o f p a r a m e t e r s f o r t h e i n t r a c e l l u l a r ODE system
		const i n t nparams = 3 9 ;
	Parameters Values Units // C e l l s G1 p h a s e a r e d e t e r m i n e d by l o w MPF c o n c e n t r a t i o n ( [MPF] < 0 . 0 6 ) Description v i // C e l l s i n M p h a s e a r e d e t e r m i n e d by : [WEE1] > [MPF] 0.025 μM/min Constant rate of cyclin synthesis const i n t index MPF = 7 ;
	v d	0.25 const i n t index WEE1 = 8 ; μM/min Maximum rate of cyclin degradation by protease X
	K d	0.025 μM // C h a r a c t e r c o n s t a n t s	Michaelis constants for cyclin degradation
	K c	0.5 // Choose t h e i n i t i a l d i s t r i b u t i o n o f c e l l s μM Michaelis constants for cyclin activation
	k d	0.01 // Type ' random ' f o r random d i s t r i b u t i o n 1/min Degradation rate of cyclin Parameters Values Units Description // Type ' f i x e d ' f o r c o n s t a n t i n i t i a l v a l u e s :
	V 1 V 2 V 3 V 4 k 1 k 2 k 3 Table B.1: Albert Goldbeter mitotic oscillator model: parameters description. 3 1/min Circadian clock // t h i s means a l l c e l l s have same i n i t i a l d i s t r i b u t i o n maximum rate of transcription 1.5 1/min c 0.01 nM Concentration of constitutive activator p 8 Unit less Hill coefficient char * i n i t i a l v a l u e s t y p e = "random" ; maximum rate of transcription 1 1/min maximum rate of transcription 0.5 1/min 0.01 unitless Michaelis constant 0.01 unitless 0.01 unitless k 7d 0.09 Degradation rate of nuclear BMAL1 h -1 k 7a 0.003 Deactivation rate of nuclear BMAL1 h -1 k 6a 0.09 Activation rate of nuclear BMAL1 h -1 k 6d 0.12 Degradation rate of nuclear BMAL1 h -1 k 6t 0.06 Nuclear export rate of BMAL1 h -1 k 5t 0.45 Nuclear import rate of BMAL1 h -1 k 5d 0.06 Degradation rate of BMAL1 h -1 k 5b 0.24 Translation rate of BMAL1 h -1 k 4d 0.75 Degradation rate h -1 k 4b 2.16 nM Michaelis constant of Bmal1 transcription k 4 0.01 unitless r 3 Unit less Hill coefficient of Bmal1 transcription Michaelis constant k 3d 0.12 h -1 Degradation rate of nuclear PER2/CRY complex ν 4b 3.6 Maximal rate of Bmal1 transcription n M -1 h -1 Michaelis constant q 2 Unit less Number of PER2/CRY2 complex forming subunits Michaelis constant k 1i 0.56 nM Inhibition constant of Per2/Cry transcription k 3t 0.02 Nuclear export rate of PER2/CRY complex h -1 k 2t 0.24 Nuclear import rate of PER2/CRY complex h -1 k 2d 0.05 Degradation rate of cytoplasmic PER2/CRY complex h -1 k 2b 0.3 n M -1 h -1 Formation rate of cytoplasmic PER2/CRY complex maximum rate of transcription ν 1b 9 Maximal rate of Per2/Cry transcription n M h -1 1 nM Michaelis constant of Per2/Cry transcription // D i v i s i o n t y p e : " d i r a c " or " k o s h l n a d " or " n o d i v i s i o n " k 1b k 1d 0.12 h -1 Degradation rate of Per2/Cry mRNA char * d i v i s i o n t y p e = " n o d i v i s i o n " ;
						Cell cycle
		k 0mpf	10	h -1		Activation rate of MPF
		k 1mpf	0.05	nM		Activation rate of MPF
	s Parameter Value Parameter Value Parameter Value Parameter Value 20 nM Inhibition constant of MPF d wee1 5 h -1 Degradation rate
	k 1 k 2 k 2 k 2 k 3	0.03 n k actw 0.03 d w1 cw 1 C k inactw 0.1 1 k 1wee1 d w2 k act	2 1 1 0.5 0 200 0.5 1 0.01	J 5 Unit less h -1 k 7 nM nM k 8 nM h -1 J 7 J 8 nM h -1 h -1	0.3 Hill coefficient k diss Activation rate of WEE1 due to BMAL1/CLOCK 0.001 V i25 1 k 13 0.1 Michaelis constant Concentration of constant activator J a25 0.25 k 14 0.1 Coupling strength to the circadian clock Deactivation rate J i25 0.001 k 15 1.5 k wee 0.001 k 16 1 k wee Michaelis constant Degradation rate Activation rate of MPF inhibitor Population model	0.25 0.01 0.01 0.15, 1.3,
	k 3 J 3	10 θ 1 θ 2 0.01 Ja J i	0.09 0.06 0.1 0.1	k 9 nM nM k 10 Unit less Unit less	0.1 Threshold value for G1 to S/G2 transition k 16 2 Threshold value for mitotic division (M to G1) 0.04 J 15 0.01 Stiffness of the switch Stiffness of the switch	k 25 k 25	0.05 5
	k 4	2		J 9		0.01	J 16	0.01	μ	0.005
	k 4	35		J 10		0.01	V awee	0.25	k 6	0.1
	J 4	0.01		k 11		0.1	V iwee	1	k 12	3
	k 5	0.005		k 12		0.01	J awee	0.01	V a25	1
	k 5	0.3		k 12		1	J iwee	0.01

  void i n i t i a l s t a t e ( double * * , double * * , const int , int , char * , double [ ] , double [ ] [ 2 ] , double [ ] , double , double ) ;

	#e n d i f
	/ * -
	/ 1D
	double * md new ( i n t ) ;
	void m d f r e e ( double * , i n t ) ;
	// 2D
	double * * md2 new ( int , i n t ) ;
	void m d 2 f r e e ( double * * , int , i n t ) ;
	i n t * * mi2 new ( int , i n t ) ;
	void m i 2 f r e e ( i n t * * , int , i n t ) ;
	double * * m d 2 c a l l o c ( int , i n t ) ;
	// 3D
	double * * * md3 new ( int , int , i n t ) ;
	void m d 3 f r e e ( double * * * , i n t , i n t , i n t ) ;
	void s t a t e s u p d a t e ( double , char * , double [ ] , const int ,
	i n t &, i n t &, i n t &, i n t &, i n t &, i n t &,
	double * * , double * * , double * * , i n t * * , i n t * * ,
	deque <int> &, deque <double> &,
	double , double , o f s t r e a m &);
	void MVfunction ( int , const int , double * * , double [ ] , deque <double> &);

// R e a l l o c a t i o n double * * m d 2 r e a l l o c ( double * * , int , int , int , i n t ) ; i n t * * m i 2 r e a l l o c ( i n t * * , int , int , int , i n t ) ; // F u n c t i o n s double f u n c t i o n G ( double , double , double , double ) ; void p r i n t p o i n t e r ( double * , const i n t ) ; void p r i n t m a t r i x ( o f s t r e a m &, double * * , const int , i n t ) ; void p a r a m s a s s i g n a t i o n ( double * , double * ) ; void p a r a m s a s s i g n a t i o n 2 ( double * , double * ) ; void O D E s y s t e m p a r a m s a s s i g n a t i o n ( double * params ) ; void r i g h t h a n d s i d e s a b i n a ( double , double * * , double * * , double * , double * * , const int , i n t ) ; void r i g h t h a n d s i d e c e l l c y c l e a l o n e ( double , double * * , double * * , double * , double * * , const int , i n t ) ; void RK4 step ( const double , double , double * * , double * * , double * , double * * , const int , i n t ) ; void d i v i s i o n t y p e ( char * , double , const double , const double , const double , int , int , double * * , i n t * * , i n t * * , int , i n t &, deque <int> &, deque <int> &, deque <double> &, deque <double> &, deque <double> &);

Acknowledgement

Appendix A Additional Work

In this Appendix, we expose additional works and results obtained during the thesis that are not directly related to the regulation of the cell cycle by the circadian clock. However, they still belong to the framework of population dynamics, structured equations and tumor.

A.1 Meningioma

A.1.1 Introduction

In this section we would like to exploit data collected for several patients who had a surgery to extract a meningioma. For each patient, we have information about the tumor size at the time of surgery, the average age of cells and the percentage of cycling cells Ki-67 index (Table A.1). Based on these data, we would like to set a model that estimates the tumor age (the time elapsed from onset until the operation date). This work was done in collaboration with Hagen Huttner, Jonas Frisén and Olaf Bergmann from the Karolinska institute in Stockholm.

To accomplish this task, we develop a mathematical model that describes the proliferation and aging of a cell population. The model consists of an age-structured partial differential equation with Gompertz growth that tracks both the tumor size and the age of cells. Gompertz law was shown to reproduce well the growth dynamics in low grade meningioma. According to it, the tumor size follows a sigmoid curve: growth is initially exponential and slows down as the tumor size approaches a maximal value [START_REF] Laird | Dynamics of tumour growth[END_REF]. During the growth, the age of tumor cells is tracked by the age-structure of the equation. We assumed that there was no cell death, that newborn cells were aged zero and that mother cells retained their age at division. Our estimates were based on two experimentally measured quantities for each

Supplementary materials

Auxiliary equations for Tyson and Novak model: -------------------------------------------------------------------p a r a m e t e r s . h ---------------------------------------------------------------------- * --------------------------------------------------o f s t r e a m w r i t e f i n a l p o p ( " f i n a l p o p . t e x t " ) ; p r i n t m a t r i x ( w r i t e f i n a l p o p , &K[ 0 ] , neqs , n n s t a t e s ) ; w r i t e f i n a l p o p . c l o s e ( ) ; //- -------------------------------------------------- -------------------------------------------------------------------- e l s e i f ( strcmp ( i n i t i a l v a l u e s , " f i x e d " ) == 0){ f o r ( i n t i = 0; i < n e q s ; i ++) { K

Appendix C

Code source files

// W r i t i n g f i n a l s t a t e //-

} // End f u n c t i o n d i v i s i o n / * ---------------------------------------------------------------------f u n c t i o n s . cpp ------------------------------------------------------------------------ *

e t e r f n c t i o n s : s w i t c h l i k e b e h a v i r : s w i t c h e s t o t h e u p p e r // s t a t e when a / b = 1 // c , d : s t i f f n e s s o f t h e s w i t c h , s t i f f e r when c , d --->

} // F u n c t i o n t h a t p r i n t s t h e v a l u e s o f a p o i n t e r

void p r i n t p o i n t e r ( double * y , const i n t n ) { f o r ( i n t i = 0; i < n ; i++ ) cout << * ( y+i ) <<e n d l ; } // F u n c t i o n t h a t p r i n t s a m a t r i x u s i n g o f s t r e a m void p r i n t m a t r i x ( o f s t r e a m &w r i t e m a t r i x , double * * K, const i n t neqs , i n t n s t a t e s ) { f o r ( i n t i = 0; i < n e q s ; i ++) { f o r ( i n t j = 0; j < n s t a t e s ; j ++) { w r i t e m a t r i x << K[ i ] [ j ] << " " ; } w r i t e m a t r i x <<" \n" ; } } // Mean v a l u e f u n c t i o n : C a l c u l a t e s t h e m o l e c u l a r a v e r a g e s o f c e l l s void MVfunction ( i n t n s t a t e s , const i n t neqs , double * * K, double params MVfunction [ ] , deque <double> &output ) { double r ; r = 0.; f o r ( i n t i = 0; i < n e q s ; i ++) { f o r ( i n t j = 0; j < n s t a t e s ; j ++) { r+=K[ i ] [ j ] ; } r * =params MVfunction [ i ] ; r = r /(double ) n s t a t e s ; output . push back ( r ) ; r = 0.; } } / * - --------------------------------------------------------------------i ------------------------------------------------------------------------ * ----------------------// i f ( strcmp ( i n i t i a t i o n t y p e , "random" ) == 0) { f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < n e q s ; i ++) { double y i 0 = i n i t s t a t e v a l u e s r [ i ] [ 0 ] ; double y i 1 = i n i t s t a t e v a l u e s r [ i ] [START_REF] Achermann | Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves[END_REF] ; K[ i ] [ j ] = r a n d o m f l o a t ( yi0 , y i 1 ) ; ----------------------// e l s e i f ( strcmp ( i n i t i a t i o n t y p e , " f i x e d " ) == 0) { f o r ( i n t j = 0; j < n s t a t e s ; j ++) { f o r ( i n t i = 0; i < n e q s ; i ++) { K[ i ] [ j ] = randomfloat ( 0 . , 1 . 1 ) ; // i n i t s t a t e v a l u e s f --------------------------------------------------------------------m d ------------------------------------------------------------------------ * -----------------------------------------------------------------------------/ -----------------------------------------------------------------------------#

n i t i a l s t a t e . cpp

//K[ i ] [ j ] = i n i t s t a t e v a l u e s f [ i ] + r a n d o m f l o a t ( 0 . , 0 . 1 ) ; } // D e f i n i n g t h e i n i t i a l c e l l c y c l e p e r i o d f o r each c e l l

/ md new : F u n c t i o n t h a t r e s e r v e s s p a c e ( d o u b l e ) f o r a v e c t o r w i t h n e l e m e n t s // mi new : F u n c t i o n t h a t r e s e r v e s s p a c e ( i n t ) f o r a v e c t o r w i t h n e l e m e n t s //-