STRUCTURAL AND PARAMETRIC IDENTIFICATION OF BAC-TERIAL REGULATORY NETWORKS

Keywords: Inférence Inference of genetic regulatory networks, algorithms for structural and parametric identification, systems biology, modelling of bacterial regulatory networks, parameter estimation vi

R ÉSUM É en français

L'interprétation des grandes quantités de données d'expression génique récemment générées par des techniques expérimentales à haut-débit exige des outils mathématiques et informatiques fiables pour l'inférence des interactions régulatrices. Nous nous intéressons à l'inférence des interactions régulatrices et l'amélioration des résultats de l'inférence en identifiant les informations précises fournies par les données expérimentales. Nous avons développé une approche expérimentale et computationnelle intégrée pour l'inférence de modèles quantitatifs de promoteurs bactériens à partir des données d'expression génique temporelles mesurée par l'intermédiaire de gènes rapporteurs fluorescents. Nous montrons comment les effets physiologiques globaux et les concentrations de protéines peuvent être estimés à partir des données de fluorescence et intégrés dans des méthodes d'inférence, à la fois structurelle et paramétrique, des fonctions de régulation génique. Nous avons validé notre approche sur un module central dans le réseau de régulation contrôlant la motilité et le système de chimiotactisme chez Escherichia coli. L'approche proposée est orthogonale aux méthodes déjà existantes pour l'inférence des réseaux de régulation à partir de données temporelles d'expression génique et peut être intégré avec plusieurs autres méthodes proposées dans la littérature.

ACKNOWLEDGEMENTS

First and foremost, I want to present my sincere gratitude to my advisor Hidde DE JONG. His full support and significant contribution during my entire PhD, his availability at any time and his kindness and patience in teaching me great science and guiding the steps of my early carrier allowed me to successfully obtain my PhD degree.

As well, thanks to Hidde, I am now emotionally attached to what is commonly only a work place, the INRIA laboratory.

My deepest acknowledgements and sympathy go to my second advisor Eugenio CINQUEMANI for his guidance, inspiring ideas and professional investment to make my PhD experience productive and stimulating and whose help had been essential for the completion of my work. I am honoured to have been his first PhD student and to learn from such a resourceful and enthusiastic teacher.

I would also like to thank my third advisor Johannes GEISELMANN for the opportunity he has given to me to learn and practice microbiology, allowing my work to become also my passion. I deeply appreciate the immense experience and pertinent advices he has brought to my thesis, but also his bright personality.

LONG ABSTRACT in English

High-throughput technologies yield large amounts of data about the steady-state levels and the dynamical changes of gene expression in bacteria. An important challenge for the biological interpretation of these data consists in deducing the topology of the underlying regulatory network as well as quantitative gene regulation functions from such data. A large number of inference methods have been proposed in the literature and have been successful in a variety of applications, although several problems remain.

We focus here on improving two aspects of the inference methods. First, transcriptome data reflect the abundance of mRNA, whereas the components that regulate are most often the proteins coded by the mRNAs. Although the concentrations of mRNA and protein correlate reasonably during steady-state growth, this correlation becomes much more tenuous in time-series data acquired during growth transitions in bacteria because of the very different half-lives of proteins and mRNA. Second, the dynamics of gene expression is not only controlled by transcription factors and other specific regulators, but also by global physiological effects that modify the activity of all genes.

For example, the concentrations of (free) RNA polymerase and the concentration of ribosomes vary strongly with growth rate. We therefore have to take into account such effects when trying to reconstruct a regulatory network from gene expression data.

We propose here a combined experimental and computational approach to address these two fundamental problems in the inference of quantitative models of the activity of bacterial promoters from time-series gene expression data. We focus on the case where the dynamics of gene expression is measured in vivo and in real time by means of fluorescent reporter genes.

Our network reconstruction approach accounts for the differences between mRNA and protein half-lives and takes into account global physiological effects. When the half-lives of the proteins are available, the measurement models used for deriving the activities of genes from fluorescence data are integrated to yield estimates of protein concentrations. The global physiological state of the cell is estimated from the activity of a phage promoter, whose expression is not controlled by any transcription factor and depends only on the activity of the transcriptional and translational machinery. We apply the approach to a central module in the regulatory network controlling motility and the chemotaxis system in Escherichia coli. This module comprises the fliA, flgM and tar genes. FliA is a sigma factor that directs RNA polymerase to operons coding for components of the flagellar assembly. The effect of FliA is counteracted by the antisigma factor FlgM, itself transcribed by FliA. The third component of the network, tar, codes for the aspartate chemoreceptor protein Tar and is directly transcribed by the FliA-containing RNA polymerase holoenzyme. The FliA-FlgM module is particularly well-suited for studying the inference problems considered here, since the network has been well-studied and protein half-lives play an important role in its functioning.

We stimulated the FliA-FlgM module in a variety of wild-type and mutant strains and different growth media. The measured transcriptional response of the genes was used to systematically test the information required for the reliable inference of the regulatory interactions and quantitative predictive models of gene regulation.

Our results show that for the reliable reconstruction of transcriptional regulatory networks in bacteria it is necessary to include global effects into the network model and explicitly deduce protein concentrations from the observed expression profiles. Our approach should be generally applicable to a large variety of network inference problems and we discuss limitations and possible extensions of the method.

R ÉSUM É SUBSTANTIEL en français

Les technologies expérimentales à haut débit produisent de grandes quantités de données sur les niveaux d'expression des gènes dans les bactéries à l'état d'équilibre ou lors des transitions de croissance. Un défi important dans l'interprétation biologique de ces données consiste à en déduire la topologie du réseau de régulation ainsi que les fonctions de régulation quantitatives des gènes. Un grand nombre de méthodes d'inférence a été proposé dans la littrature. Ces méthodes ont été utilisées avec succès dans une variété d'applications, bien que plusieurs problèmes persistent.

Nous nous intéressons ici à l'amélioration de deux aspects des méthodes d'inférence.

Premièrement, les données transcriptomiques reflètent l'abondance de l'ARNm, tandis que, le plus souvent, les composants régulateurs sont les protéines codées par les ARNm. Bien que les concentrations de l'ARNm et de protéines soient raisonnablement corrélées à l'état stationnaire, cette corrélation devient beaucoup moins évidente dans les données temporelles acquises lors des transitions de croissance à cause des demi-vies très différentes des protéines et des ARNm. Deuxièmement, la dynamique de l'expression génique n'est pas uniquement contrôlée par des facteurs de transcription et d'autres régulateurs spécifiques, mais aussi par des effets physiologiques globaux qui modifient l'activité de tous les gènes. Par exemple, les concentrations de l'ARN polymérase (libre) et les concentrations des ribosomes (libres) varient fortement avec le taux de croissance. Nous devons donc tenir compte de ces effets lors de la reconstruction d'un réseau de régulation à partir de données d'expression génique.

Nous proposons ici une approche expérimentale et computationnelle combinée pour répondre à ces deux problèmes fondamentaux dans l'inférence de modèles quantitatifs de promoteurs bactériens à partir des données temporelles d'expression génique. Nous nous intéressons au cas où la dynamique de l'expression génique est mesurée in vivo et en temps réel par l'intermédiaire de gènes rapporteurs fluorescents. Notre approche d'inférence de réseaux de régulation tient compte des différences de demi-vie entre l'ARNm et les protéines et prend en compte les effets physiologiques globaux. Lorsque les demi-vies des protéines sont connues, les modèles expérimentaux utilisés pour dériver les activités des gènes à partir de données de fluorescence sont intégrés pour estimer les concentrations des protéines. L'état physiologique global de la cellule est estimé à partir de l'activité d'un promoteur de phage, dont l'expression n'est contrôlée par aucun des facteurs de transcription et ne dépend que de l'activité de la machinérie dexpression génique.

Nous appliquons l'approche à un module central dans le réseau de régulation contrôlant la motilité et le système de chimiotactisme chez Escherichia coli. Ce module est compose des genes fliA, flgM et tar. FliA est un facteur sigma qui dirige l'ARN polymérase vers les opérons codant pour des composants de l'assemblage des flagelles. L'effet de FliA est contrecarré par le facteur anti-sigma FlgM, lui-même transcrit également par FliA. Le troisième composant du réseau, tar, code pour la protéine récepteur chimiotactique de l'aspartate, Tar, et est directement transcrit par FliA associé à l' holoenzyme ARN polymérase. Le module FliA-FlgM est particulièrement bien adapté pour l'étude des problèmes d'inférence considérés ici, puisque le réseau a été bien étudié et les demivies des protéines jouent un rôle important dans son fonctionnement.

Nous avons stimulé le module FliA-FlgM dans une variété de souches de type sauvage et mutantes et dans des milieux de croissance différents. La réponse transcriptionnelle des gènes mesurée a été utilisée pour tester systématiquement les informations requises pour l'inférence fiable des interactions régulatrices et des modèles prédictifs quantitatifs de la régulation des gènes.

Nos résultats montrent que, pour la reconstruction fiable de réseaux de régulation transcriptionelle chez les bactéries, il est nécessaire d'inclure les effets globaux dans le modèle de réseau et d'en déduire de manière explicite les concentrations des protéines à partir des profils d'expression observés, car la demi-vie de l'ARNm et des protéines sont très différentes. Notre approche reste généralement applicable à une grande variété de problèmes d'inférence de réseaux et nous discutons les limites et les extensions possibles de la méthode.

x

Bacteria in their environment

Bacteria are single-cell organisms that, under favorable conditions such as nutrient abundance, can grow and divide rapidly. The number of bacteria on earth is estimated to be 4 -6 • 10 30 [START_REF] Whitman | Prokaryotes: the unseen majority[END_REF] and their total amount of carbon or biomass equals the estimated total carbon in plants. Typically few micrometers in dimension, these numbers are indicative of their pervasiveness in our world. Their contribution is everywhere enormous, in soils, water and air, recycling all types of matter. Bacteria not only inhabit human organisms, we even rely on them to metabolize and absorb essential nutrients, to fight pathogens, and to train and improve our immune function [START_REF] Peterson | The NIH Human Microbiome Project[END_REF][START_REF] Turnbaugh | The human microbiome project[END_REF].

Bacteria can adapt to many different environments, for example they survive extreme stress conditions associated with high concentrations of toxic metals or radioactive environments [START_REF] Daly | Accumulation of Mn(II) in Deinococcus radiodurans facilitates gamma-radiation resistance[END_REF][START_REF] Keller | Tapping into microbial diversity[END_REF][START_REF] Schaechter | [END_REF].

Recent studies of the human microbiome have begun to characterize the bacterial communities living in our gastrointestinal tract [START_REF] Peterson | The NIH Human Microbiome Project[END_REF]. These communities thrive in the extreme conditions they have to deal with.

One of the microbes of the human gut flora, which can be grown easily and inexpensively in a laboratory setting, is Escherichia coli. The bacterium produces menaquinones [START_REF] Conly | The contribution of vitamin K2 (menaquinones) produced by the intestinal microflora to human nutritional requirements for vitamin K[END_REF][START_REF] Suvarna | Menaquinone (vitamin K2) biosynthesis: localization and characterization of the menA gene from Escherichia coli[END_REF], known as K2 vitamin, which plays an important and complex role in hemostasis. However, some serotypes have developed ways to thwart the immune system and can cause serious food poisoning in their hosts [START_REF] Brown | Revisiting the host as a growth medium[END_REF][START_REF] Tenaillon | The population genetics of commensal Escherichia coli[END_REF]. Escherichia coli has served for over 60 years as a model organisms for microbiology and biotechnology research [START_REF] Keseler | EcoCyc: fusing model organism databases with systems biology[END_REF][START_REF] Salgado | RegulonDB v8.0: omics data sets, evolutionary conservation, regulatory phrases, crossvalidated gold standards and more[END_REF]. Recently, this bacterium has been engineered to create "living 1. INTRODUCTION materials", assembling rows of gold nanowires. The resulting network conducts electricity and could be worth exploring for use in energy applications, such as batteries and self-healing materials [START_REF] Chen | Synthesis and patterning of tunable multiscale materials with engineered cells[END_REF].

Responses of bacteria to changes in their environment

Bacteria have developed many original solutions to respond to often rapid changes in their environment. Bacterial stress responses enable cells to survive adverse and fluctuating conditions in their environment, such as the depletion of nutrients, changes in pH and temperature, high population density [START_REF] Storz | Bacterial stress responses[END_REF]. When nutrients become exhausted, cell membranes become thicker for protection and cell division is dramatically reduced or turned off to prevent energy expenditure.

Another example of a bacterial stress response is motility [START_REF] Berg | coli in motion[END_REF]. This is thought to be one of the most impressive evolutionary aspects of bacterial behavior, as it allows bacteria to populate regions rich in nutrients and to avoid repellents [START_REF] Berg | coli in motion[END_REF]. Many bacterial species swim by rotating external filamentous organelle, known as flagella. The flagellum has three main components (as described by [START_REF] Berg | coli in motion[END_REF]), a basal body integrated in the cell wall (''the motor''), a short joint structure (''the hook'') and a long filament (''the propeller''). In E. coli, each flagellum is driven at its base by a reversible flagellar motor that is powered by a chemical gradient across the membrane and propels bacteria in a particular direction. When flagellar rotation is counterclockwise, the flagella push the bacterium forward, allowing a reasonably smooth "run". When the rotation is clockwise, the flagella pull in opposing directions and the bacteria "tumble" [START_REF] Berg | coli in motion[END_REF] and Figure 1.1). Sensing chemicals is achieved through a complicated protein chemotaxis system [START_REF] Porter | Signal processing in complex chemotaxis pathways[END_REF] that controls the flagellar motor, allowing bacteria to migrate towards environments that are optimal for growth.

If the movement is directed towards an attractant, the running period is prolonged and bacteria swim towards the attractant. On the contrary, if a repellent is encountered, the bacteria tumble, which prevents them from swimming towards the repellent.

Context

Figure 1.1: Bacterial motility. Most of bacterial species can move by rotating their flagella when sensing chemical gradients. The complex system coordinating movement and chemical sensing is the chemotaxis system. Cells switch continuously between two modes of movement: a "run" mode in which the bacterium swims forward, by rotating their flagella counterclockwise and a "tumble" mode in which the bacterium randomly changes direction, by rotating their flagella clockwise. The rotary motor located at the base of each flagellum controls the direction of rotation [START_REF] Berg | coli in motion[END_REF]. The complex chemotaxis system, triggered by signals emitted by the transmembrane receptor proteins, dictates the direction of rotation of the flagella.

Responses controlled by complex regulatory networks

How does a bacterial cell initiate and coordinate its adaptive responses when changes in their environment occur? This is achieved by highly sophisticated, complex regulatory networks.

For example, the synthesis and function of the flagellar and motility system is based on the coordinated expression of more than 50 genes [START_REF] Chilcott | Coupling of flagellar gene expression to flagellar assembly in Salmonella enterica serovar typhimurium and Escherichia coli[END_REF]. The expression of these genes responds to environmental stimuli and, in addition, to signals that are coupled to the morphological development of the flagella. The flagellar genes are organized in a transcriptional hierarchy of three operon classes [START_REF] Chevance | Coordinating assembly of a bacterial macromolecular machine[END_REF][START_REF] Kutsukake | Transcriptional analysis of the flagellar regulon of Salmonella typhimurium[END_REF]Macnab, 1996a) as shown in Figure 1.2. The class 1 operon, flhDC, encodes the proteins FlhC and FlhD, which form a heteromultimeric complex initiating the transcription of the entire flagellar cascade through the class 2 operons. These operons encode the structural proteins required for flagellar

INTRODUCTION

hook assembly as well as the main regulator of the class 3 operons, the sigma factor FliA (σ 28 ). When bound to core RNA polymerase, FliA directs the transcription of the class 3 operons, the lower level of the regulatory hierarchy. However, FliA activity is inhibited by the anti-sigma factor FlgM, which binds to FliA and thus prevents its association with RNA polymerase, delaying class 3 genes expression and completion of the flagellar assembly. FlgM is transcribed from both a class 2 and class 3 promoter and can be excreted once the intermediary hook basal body structure is constructed.

Cellular processes are very complex, but it seems that such processes can often be broken down into a limited number of recurring patterns of connectivity. The transcriptional cascade in Figure 1.2 can be decomposed into elementary network motifs, such as the SUM input FeedForwardLoop [START_REF] Alon | Network motifs: theory and experimental approaches[END_REF]. That is, the master regulator FlhDC activates a second regulator, FliA, and both activate, in an additive fashion, the operons that produce the flagella motor. This specific network motif prolongs flagella expression following deactivation of the master regulator, highlighting the regulatory hierarchy and timing of the control of the flagella assembly [START_REF] Kalir | Using a quantitative blueprint to reprogram the dynamics of the flagella gene network[END_REF][START_REF] Kalir | A coherent feed-forward loop with a SUM input function prolongs flagella expression in Escherichia coli[END_REF].

Experimental measurements

Recent advances in molecular biology and in biophysics have led to new technologies for measuring cellular processes at the molecular level and in real time, including DNA microarrays and RNA sequencing, gene reporter systems, quantitative RT-PCR and mass-spectrometry based measurement of proteins. This allows the stress responses of bacteria to be monitored and may give insights into the functioning of the regulatory networks controlling these responses.

Measurements of the transcriptome of a bacterial cell by means of DNA microarrays

or RNA sequencing produce quantitative information on the state of the entire transcriptional program of an organism at any time during an experiment [START_REF] Dharmadi | Dna microarrays: experimental issues, data analysis, and application to bacterial systems[END_REF]. These approaches provide a relative quantification of mRNA abundance and cannot be obtained in vivo, in the sense that the actual measurements are carried out on molecules extracted from the cells. Techniques such as quantitative RT-PCR also measure relative mRNA concentrations, though usually not on a genome-wide scale [START_REF] Saunders | Real-Time PCR: Advanced Technologies and Applications[END_REF], see [START_REF] White | High-throughput microfluidic single-cell RT-qPCR[END_REF] for an exception). Fluorescent The hierarchical transcriptional network controlling flagellar assembly and the synthesis of the chemotaxis system (Karlinsey et al., 2000a). Class 1 genes encode the proteins FlhD and FlhC, transcriptional activators that are required for expression of all the other flagellar genes. Class 2 genes encode structural components of the Hook-Basal Body (HBB) complex, as well as regulatory proteins, including the sigma factor σ 28 and the anti-σ 28 factor FlgM. σ 28 is required for the expression of all class 3 gene promoters, including those for flagellin and those related to chemotaxis and motility. FlgM binds to and inhibits σ 28 until completion of the HBB. FlgM is then secreted from the cell, and σ 28 -dependent (class 3) gene expression is initiated.

reporter genes allow measuring the activities of the promoter regulating the transcription of genes in vivo and in real time. Although quantitative proteomics has much advanced recently [START_REF] Picotti | Selected reaction monitoring-based proteomics: workflows, potential, pitfalls and future directions[END_REF], in vivo and real time measurements of proteins on a genome-wide scale are not yet possible.

Reporter gene technology, in particular, allows the direct observation/quantification of

INTRODUCTION

the activity of genes in the cell. For motility gene data obtained by means of fluorescent reporter systems this reflects the expression of the different classes of genes in the flagellar cascade [START_REF] Kalir | Using a quantitative blueprint to reprogram the dynamics of the flagella gene network[END_REF]. We can observe the precise temporal expression profiles of the most important genes involved in flagellar assembly. However, discovering the structure and functioning of the regulatory network is not directly possible from the raw experimental data.

Network inference

Computational tools are needed to uncover regulatory interactions from the large amounts of experimental observations, as well as to construct dynamic models of the functioning of bacteria. This problem is commonly defined as reverse engineering [START_REF] Tegner | Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling[END_REF], network reconstruction [START_REF] Maccarthy | Using large-scale perturbations in gene network reconstruction[END_REF] or network inference [START_REF] Faith | Reverse-engineering transcription control networks[END_REF][START_REF] Gardner | Inferring genetic networks and identifying compound mode of action via expression profiling[END_REF].

Different modeling formalisms have been proposed for inferring the topology of gene regulatory networks from gene expression data ( [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF][START_REF] Hecker | Gene regulatory network inference: data integration in dynamic models-a review[END_REF][START_REF] Villaverde | Reverse Engineering Cellular Networks with Information Theoretic Methods[END_REF]. These models represent regulatory networks in different ways, for example, as (oriented) graphs (Bayesian networks), discrete dynamic systems (Boolean networks) and continuous dynamic systems (differential equations). Specific inference algorithms for each of these formalisms have been developed, reconstructing the activity of one gene as a function of the activity of other genes. Although this has resulted in powerful methods, a number of recurrent problems remain [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF][START_REF] Prill | Towards a rigorous assessment of systems biology models: the dream3 challenges[END_REF].

Examples of such issues are the dimensionality problem, or the problem of discriminating direct from indirect regulations. Dimensionality issues arise, e.g. when the number of variables (genes) is much larger than the number of experimentally observed quantities. This is typical for in vivo time course gene expression measurements for quantitative dynamic model inference and results in a multitude of regulatory structures consistent with the data. Moreover, measurements are typically noisy and do not capture the entire dynamics of the gene expression. Using additional constraints, biological background knowledge and clever simplifications can help analyze the prohibitively large space of possible solutions in a time-efficient manner (de [START_REF] Smet | Advantages and limitations of current network inference methods[END_REF]. Second, a regulator, e.g., a transcription factor may directly control 1.2 Problem Statement: two problems in network inference the expression of a target gene by fixing to its promoter region, but can also affect it indirectly, by regulating the expression of an intermediary gene whose product in turn regulates the expression of the target gene. Such direct and indirect regulatory influences are not easily distinguishable from gene expression data [START_REF] Basso | Reverse engineering of regulatory networks in human B cells[END_REF][START_REF] Rice | Reconstructing biological networks using conditional correlation analysis[END_REF]. This is thought to be due to the fact that the inference algorithms rely on specific assumptions about the underlying network topology (e.g. cyclic or acyclic network structures, including feedback loops or cascades). However, network motif analysis makes it possible to quantitatively assess how the difficulty of distinguishing direct and indirect connections affects inference methods [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF]. These different problems may lead to systematic errors in predicting regulatory interactions and, as a result, compromise the performance of network inference algorithms [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF].

Several other fundamental problems of this sort exist in model inference. In my PhD work I have focused on two other recurrent problems in network inference that will be introduced in the next section.

Problem Statement: two problems in network inference

Usually, network inference algorithms rely on transcriptome data generated for instance, through DNA microarray analysis or RNA sequencing. These relative RNA concentrations characterize the transcriptome state of the cell as well as the activity of the promoters initiating the transcription of the genes. The problem of inferring regulatory interactions is that in general the active regulator is not mRNA but protein.

At steady-state, mRNA and protein concentrations are relatively well correlated [START_REF] Lu | Absolute protein expression profiling estimates the relative contributions of transcriptional and translational regulation[END_REF][START_REF] Taniguchi | Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells[END_REF]. However, this is not expected to occur when the two vary over time. mRNA and protein have different half-lives and their concentrations evolve on different time-scales. For instance, in bacteria, mRNA half-lives are on the order of a few minutes [START_REF] Bernstein | Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays[END_REF], whereas most of the proteins are stable [START_REF] Larrabee | The relative rates of protein synthesis and degradation in a growing culture of Escherichia coli[END_REF][START_REF] Mosteller | Metabolism of individual proteins in exponentially growing Escherichia coli[END_REF] and the degradation rate is dominated by growth dilution, taking place on the time-scale of a cell cycle. The effect of rapid responses in gene expression, taking place within a single generation, may thus persist

INTRODUCTION

over several generations [START_REF] Maier | Correlation of mRNA and protein in complex biological samples[END_REF][START_REF] Taniguchi | Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells[END_REF]. As a consequence of this temporal decorrelation of mRNA and protein concentrations, inference of regulatory networks relying exclusively on time-series transcriptome data may potentially lead to spurious results.

Another important problem derives from the fact that the dynamics of gene expression is not only controlled by transcription factors, small regulatory RNAs, and other specific regulators, but also by global physiological effects influencing the rates of transcription and translation of all genes (Berthoumieux et al., 2013b;[START_REF] Gerosa | Dissecting specific and global transcriptional regulation of bacterial gene expression[END_REF][START_REF] Keren | Promoters maintain their relative activity levels under different growth conditions[END_REF][START_REF] Klipp | Growth-rate-dependent partitioning of RNA polymerases in bacteria[END_REF]. Most gene expression studies have been based on the assumption that cells produce similar levels of total RNA per cell, without including standardized controls that would reveal global transcriptional amplification or repression. For instance, cells can globally up-regulate their gene expression program, producing two to three times more total RNA and generating larger cells. In the conventional approach to expression analysis, normalized amounts of RNA would be introduced into the assay, thus masking changes in the activity of gene expression machinery. Ignoring such changes, for example in experiments with important variations of the growth rate, can lead to artefacts in inferring regulatory interactions [START_REF] Lovén | Revisiting global gene expression analysis[END_REF][START_REF] Regenberg | Growth-rate regulated genes have profound impact on interpretation of transcriptome profiling in Saccharomyces cerevisiae[END_REF]. Unfortunately, global physiological parameters characterizing changes in the activity of the gene expression machinery, such as the concentrations of (active) RNA polymerase and ribosome, are difficult to quantify in a direct way.

The aim of my PhD work is (I) to propose a combined experimental and computational approach to address the above two fundamental problems in the inference of quantitative models of regulatory interactions from time-series data and (II) the application of these methods to real data of gene expression from the regulatory network controlling motility in E. coli. This network has been well-studied and is therefore particularly suitable as a test-case. From a methodological point of view, several issues will be addressed during the project, such as the choice of appropriate modeling formalisms, the study and integration of time-series data in dynamic models, the design of informative experiments, the application of effective algorithms for model identification on real data, and the interpretation of the identification results to learn about the structure and functioning of the network. In order to generate a set of rich observations [START_REF] Giepmans | The fluorescent toolbox for assessing protein location and function[END_REF][START_REF] Southward | The dynamic microbe: Green fluorescent protein brings bacteria to life[END_REF], since they allow gene expression to be monitored with high precision and temporal resolution. Fluorescent reporter genes consist of transcriptional fusions of a gene encoding a fluorescent protein, e.g., GFP or mCherry, to the promoters of the target genes, on (low-copy) plasmids or on the chromosome. Chromosomal reporters avoid a number of potential artifacts, such as a change in plasmid copy number across different growth phases, but they are more difficult to construct and the intensity of the fluorescence signal may be close to the background fluorescence, especially when GFP is used. In this study, we have chosen to use plasmidic reporters of the motility genes, available in a reporter library [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF]. The copy number of these plasmids was shown previously to be stable in media supporting different growth rates (Berthoumieux et al., 2013b;[START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF].

Reconstruction of protein concentrations and global regulatory effects

The experimental protocol allows monitoring the transcriptional response of a biological network. How can we reconstruct quantities of interest for our purpose from these data, notably protein concentrations and global physiological effects?

In the case of fluorescent reporter gene system, primary absorbance and fluorescence signals can be transformed into what are commonly called in the literature promoter activities [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF], using kinetic measurement models of gene expression (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF][START_REF] Huang | Integrated modeling and experimental approach for determining transcription factor profiles from fluorescent reporter data[END_REF][START_REF] Ronen | Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics[END_REF][START_REF] Wang | Mathematical analysis and quantification of fluorescent proteins as transcriptional reporters[END_REF]. More

INTRODUCTION

precisely, reporter gene data allow one to deduce protein synthesis rates and under certain conditions to consider them proportional to mRNA concentrations and promoter activities. These quantities reflect the transcriptional activity of the gene, but they can also be used to reconstruct the concentration of the protein, using information on the protein half-life (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF].

As explained in Section 1.2, large-scale differences in gene expression over time or across conditions may also reflect global changes in cellular physiology. The approach described in Berthoumieux et al. (2013b) allows the global state of the cell to be monitored in real-time and in vivo during the growth transition. I will therefore use a GFP reporter driven by a constitutive promoter, not regulated by any transcription factor, to assay the time-varying physiological state of the cell. For example, a plasmid expressing a GFP reporter for a phage promoter [START_REF] Oppenheim | Switches in bacteriophage lambda development[END_REF], not regulated by any protein in the host cell, can provide this type of information. The variations in the activity of the constitutive promoter reflect changes in the overall physiological state of the cell, including the RNA polymerase and ribosome concentrations, as well as pool sizes of amino acids and nucleotides.

Transcriptional response of the FliA-FlgM module

We apply the experimental approach to a central module in the regulatory network controlling the synthesis of flagella and the chemotaxis sensing system in Escherichia coli [START_REF] Chevance | Coordinating assembly of a bacterial macromolecular machine[END_REF][START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF]Macnab, 1996a). This module comprises the FliA and FlgM transcription factors and their targets. FliA or σ 28 is a sigma factor that directs RNA polymerase to operons coding for the flagellar filament and the chemotaxis sensing system controlling the flagellar motor. The effect of FliA is counteracted by the anti-sigma factor FlgM. As typical examples of FliA-dependent genes we study flgM, the gene encoding FlgM, and tar. The latter gene encodes the aspartate chemoreceptor protein Tar, which activates the flagellar motor component [START_REF] Berry | Microbiology. How bacteria change gear[END_REF][START_REF] Macnab | Escherichia coli and Salmonella: Cellular and Molecular Biology, chapter Chemotaxis[END_REF]. The FliA-FlgM module forms a checkpoint in the temporally-organized expression cascade. It is particularly well-suited for investigating the inference problems considered here, since the interactions in this network have been well-studied and protein stability has been found to play an important role in its functioning.

Principal research questions and approaches

We will experimentally excite the FliA-FlgM module in a variety of wild-type and mutant conditions, in different growth media, and measure the transcriptional response of the genes. Promoter activities and protein concentrations will be computed for the fliA, flgM and tar genes, as well as the activity of the constitutive phage promoter pRM, to account for global regulatory effects.

Inference of quantitative models of the FliA-FlgM module

How can the data on the transcriptional response of the FliA-FlgM module be used to systematically test the information required for the reliable inference of the regulatory interactions (structure) and quantitatively predictive models (parameters) of gene regulation?

For the structural inference problem I will use a previously described inference method (Porreca et al., 2010a). I will notably test if the use of fliA and flgM promoter activities, instead of their protein concentrations, allows us to retrieve the expected pattern of regulatory interactions. Furthermore, I will assess to which extent the results can be improved when considering protein concentrations instead of promoter activities.

In addition, analysis will consider the presence or absence of a model factor for the potential effects of global physiology. In order to compute protein concentrations we will use measured or estimated half-lives, while global physiological effects will be measured by means of a constitutive promoter.

In order to quantify to which extent the addition of the latter information improves the identification of quantitative models of promoter activity, I will construct kinetic models of the regulation of FliA-dependent genes. Using heuristic global optimization methods, such as genetic algorithms, I will then estimate parameter values from the data and asses the quantitative fit in the different conditions.

INTRODUCTION

Contributions

The interpretation of the large amounts of data yielded recently by high-throughput experimental techniques requires more reliable mathematical and computational tools for the inference of regulatory interactions.

In this work, we have made explicit the relation between experimental data and physiological quantities by means of mathematical models of gene expression, calling into question two assumptions that are commonly made in the inference of regulatory interactions and quantitative gene regulation functions from time-series data. The first assumption is that transcriptome data alone are sufficient to capture the time-varying state of gene expression. Often, the regulators of gene expression are proteins, while mRNA and protein concentrations are not correlated in dynamic experiments. As a consequence, currently it is not possible to fully exploit the information contained in time-series transcriptome data [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF]. A second implicit assumption in the analysis of transcriptome data is that gene regulation can be reduced to the action of transcription factors and other specific regulators. This ignores the fact that the activity of the transcriptional and translational machinery, as well as other global physiological effects, may drastically change over the course of an experiment, a fact that has been well-documented for microorganisms [START_REF] Dennis | Control of rRNA synthesis in Escherichia coli: a systems biology approach[END_REF][START_REF] Scott | Bacterial growth laws and their applications[END_REF]. This may lead to erroneous interpretations and the inference of spurious regulatory interactions [START_REF] Lovén | Revisiting global gene expression analysis[END_REF].

The main contribution of this thesis is an integrated experimental and computational approach for addressing the above two problems, in the context of time-series measurements of gene expression by means of fluorescent reporter genes. We notably show how global physiological effects and protein concentrations can be estimated from fluorescence data and integrated into methods for the inference of structural and parametric gene regulation functions. This work relies on solid results obtained previously. The reconstruction of protein concentrations from real-time promoter activities by means of kinetic models as well as the quantification of global physiological effects by means of reporter genes have been proposed before (Berthoumieux et al., 2013b;[START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF][START_REF] Gerosa | Dissecting specific and global transcriptional regulation of bacterial gene expression[END_REF][START_REF] Keren | Promoters maintain their relative activity levels under different growth conditions[END_REF]. For instance, [START_REF] Gerosa | Dissecting specific and global transcriptional regulation of bacterial gene expression[END_REF] have developed quantitative models to dissect global and specific regulation of E. coli genes involved in arginine biosynthesis [START_REF] Gerosa | Dissecting specific and global transcriptional regulation of bacterial gene expression[END_REF].

Contributions

To our knowledge, the work presented here is the first systematic study of how the integration of information on both global physiological effects and protein concentrations can improve the inference of regulatory interactions and the identification of gene regulation functions from time-series data. It is important to emphasize that the proposed approach is orthogonal to existing methods for the inference of regulatory networks from time-series gene expression data and can be combined with any of the large variety of methods proposed in the literature.

We have validated our approach by analyzing a central module of the motility network in E. coli. The FliA-FlgM module has been very well-studied and has characteristics that make it particularly suitable for our purpose, such as short half-lives due to export of certain proteins from the cell and proteolysis. The secretion and degradation rates change across conditions, depending on the strength of induction of the flagella synthesis network. As a consequence, the FliA and FlgM concentrations are expected to vary during the course of an experiment and across the experimental conditions. This yields a rich and challenging data set for testing how accounting for the distinction between cellular responses on the level of mRNA and protein influences the results of the inference process.

Furthermore, we use the reporter gene data not only for deducing the regulatory structure but also for quantifying the regulation function of two FliA-dependent motility genes, not known to be regulated by any other transcription factors. When progressively solving the problems mentioned above, by integrating information about the activity of the gene expression machinery and computing estimates of protein concentrations from promoter activities, both the structure and the dynamics of the regulation of the tar and flgM promoters could be identified successfully. We emphasize that, when using available measurements of FliA and FlgM half-lives, this was achieved without increasing the number of parameters in the models and is therefore not simply a consequence of increasing the degrees of freedom. The results confirmed the important roles played by global physiological effects and the active regulation of FliA and FlgM half-lives in shaping the dynamics of FliA-dependent promoters.

We believe that the approach proposed in this work has broad practical applicability for exploiting and analyzing transcriptome data and improving network inference in a variety of organisms.

1. INTRODUCTION

Thesis overview

The manuscript of this thesis will be organized as follows:

• Chapter 2 (State of the art) will review existing inference methods for dynamical models. We will also describe the experimental techniques that can be currently used to monitor gene expression data in microbial cells. Finally, data analysis methods that allow the transformation of reporter gene data into biologically relevant quantities will be presented.

• Chapter 3 (Results) will present the results obtained during this thesis. We will investigate how reconstructing protein concentrations from promoter activities and monitoring the global physiological state of the cell may improve the structural and parametric inference of gene regulatory networks. The approach will be exemplified by means of a central module of the motility network in E. coli .

A paper containing the work presented in this chapter is in preparation.

• Chapter 4 will summarize the conclusions drawn from the current work and will present perspectives and future improvements for the inference of quantitative models of regulatory networks in biology.

State of the Art

Inference of gene regulatory networks generally deals with the problem of reconstructing interactions among genes from experimental data. Different types of data exist, and correspondingly, different interaction models and different methods for their inference are most adapted to the data considered. Depending on the context, the problem goes under alternative names, notably reverse-engineering and network identification.

In this chapter we present the various methods employed in the literature for inferring gene regulatory networks. In Section 1 we introduce common experimental techniques that provide high-throughput measurements of regulatory molecules such as mRNAs or proteins. Then, we present the inference algorithms that have been developed to exploit these experimental data, their strengths and weaknesses. In the last section (Section 3), we address the data processing methods and the development of measurement models for reporter gene techniques.

Experimental techniques for measuring gene expression

Measurements of the transcriptome and proteome of bacterial cells by means of DNA microarrays, RNA sequencing, and other high-throughput or quantitative technologies have created huge amounts of data on the state of the transcriptional program in different growth conditions and genetic backgrounds, over the time course of an experiment.

Large efforts have been made to develop such experimental methods. This section offers a short overview of the experimental technologies as well as examples of their application.

High-throughput transcriptomics

Transcriptomics allows the monitoring of the genome-wide transcriptional response of the cell to an environmental stimulus or a genetic modification.

STATE OF THE ART

The most common technique used for producing data for the inference of gene networks are microarrays. There are several types of microarrays, but DNA microarrays are by far the most widely used.

The principle of DNA microarrays is based on the complementarity of messenger RNA (mRNA) and the DNA strand from which it is transcribed. Therefore mRNA will bind to single -stranded DNA molecules with the same sequence as the originally transcribed gene.

To determine gene expression patterns in a cell, the messenger RNA molecules in a sample are extracted. Each mRNA molecule is then reverse transcribed by a reverse transcriptase (RT) thus generating a complementary DNA (cDNA). In addition, the cDNA is labeled in the process, typically using a fluorescent nucleotide. Next, the labeled cDNAs obtained from cellular target mRNA are added to the microarray where they hybridize with their complementary single-stranded DNA probe fixed on the microarray substrate [START_REF] Dharmadi | Dna microarrays: experimental issues, data analysis, and application to bacterial systems[END_REF]. The fluorescence intensity measured at a particular spot in the array reveals the amount of the gene transcript which was present in the sample. Many choices of DNA microarray platforms are available, such as cDNA microarrays or oligonucleotide microarrays [START_REF] Lockhart | Expression monitoring by hybridization to high-density oligonucleotide arrays[END_REF][START_REF] Marshall | Getting the noise out of gene arrays[END_REF][START_REF] Schena | Quantitative monitoring of gene expression patterns with a complementary DNA microarray[END_REF]. Starting with the dot-blots [START_REF] Southern | Molecular interactions on microarrays[END_REF], DNA microarrays have evolved to filter arrays (nylon membrane support), and to a glass slides format, which has the advantage of high-probe density [START_REF] Dharmadi | Dna microarrays: experimental issues, data analysis, and application to bacterial systems[END_REF].

Informative high-throughput datasets have been obtained by means of DNA microarrays. Many studies on bacterial genetics make use of these datasets, especially for characterizing bacterial responses to environmental changes, trying to explore the highly complex regulatory networks and transcriptional regulation or genetic and metabolic engineering (see review of [START_REF] Dharmadi | Dna microarrays: experimental issues, data analysis, and application to bacterial systems[END_REF]). For instance, global transcriptional profiling of E. coli in acetate cultures was investigated using DNA microarrays on glass slides [START_REF] Oh | Global expression profiling of acetate-grown Escherichia coli[END_REF]. Other examples of the use of DNA microarrays datasets for E. coli can be found in the literature, such as the identification of regulatory networks [START_REF] Faith | Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF] or metabolic engineering [START_REF] Park | Global physiological understanding and metabolic engineering of microorganisms based on omics studies[END_REF].

A limitation of DNA microarray technology is that it only measures the relative concentration of mRNA. Moreover, while microarrays have contributed to our understanding

Experimental techniques for measuring gene expression

of transcription regulation (van Vliet, 2010), they have a limited precision of measurements due to variability from the array fabrication process, systematic errors in the hybridization process and heterogeneity of experimental design procedures [START_REF] Dharmadi | Dna microarrays: experimental issues, data analysis, and application to bacterial systems[END_REF]. This also increases the difficulty of combining and interpreting different available microarray datasets [START_REF] Bloom | Measuring differential gene expression by short read sequencing: quantitative comparison to 2-channel gene expression microarrays[END_REF].

While DNA microarrays thus have limitations in their applicability, sequencing technologies recently became available for the detection and quantification of transcripts in microorganisms. In addition to measuring mRNA levels, RNA-seq technology has also been used for identifying small regulatory RNAs [START_REF] Waters | Regulatory RNAs in bacteria[END_REF]. For example, [START_REF] Perkins | A strand-specific RNA-Seq analysis of the transcriptome of the typhoid bacillus Salmonella typhi[END_REF] in a study on Salmonella enterica serovar Typhi, used RNA-seq information to identify novel noncoding RNA sequences and new members of regulons. A limitation of current sequencing technologies is that the mRNA has to be extracted from a bacterial sample and measurements cannot be made in vivo or used to investigate single-cell dynamics. If one wants to evaluate a dynamic response, sampling time will become critical for the identification of changes in gene expression.

Quantitative proteomics

DNA microarrays are used to estimate mRNA levels, a molecule that rapidly changes its concentration in response to regulatory signals. However, the biologically active regulator is usually not mRNA, but protein [START_REF] Cox | Quantitative, high-resolution proteomics for data-driven systems biology[END_REF]. Although in steady-state condition proteins and mRNA concentrations are moderately correlated [START_REF] Lu | Absolute protein expression profiling estimates the relative contributions of transcriptional and translational regulation[END_REF][START_REF] Taniguchi | Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells[END_REF], if one wants to evaluate changes in gene expression dynamically over time, this is not expected to be the case. Mass spectrometry (MS)-based proteomics has emerged as an universal method for the measurement of proteins [START_REF] Bensimon | Mass spectrometry-based proteomics and network biology[END_REF][START_REF] Schmidt | Directed mass spectrometry: towards hypothesis-driven proteomics[END_REF][START_REF] Wepf | Quantitative interaction proteomics using mass spectrometry[END_REF]. Protein samples are extracted from cells and digested into peptides. The resulting peptide mixture is separated, typically by high-performance liquid chromatography (HPLC) and converted to gas phase ions by using the electrospray [START_REF] Fenn | Electrospray ionization for mass spectrometry of large biomolecules[END_REF] or matrix-assisted laser ionization [START_REF] Hillenkamp | Matrix-assisted laser desorption/ionization mass spectrometry of biopolymers[END_REF] methods. Next, the mass spectrometer scans the entire mass range every few seconds. The data analysis software then isolates the selected peptides in the mass spectra, fragments them and measures the mass spectra of each of the fragments with a high resolution. Peptide-based proteomics does

STATE OF THE ART

not directly identify proteins, but reconstructs them from the obtained mass spectra of the fragments [START_REF] Nesvizhskii | Interpretation of shotgun proteomic data: the protein inference problem[END_REF]. Using this methodology, the first complete proteome quantification for yeast has recently been published [START_REF] De Godoy | Comprehensive massspectrometry-based proteome quantification of haploid versus diploid yeast[END_REF]. One of the main perspectives is to combine proteomics and data provided by other high-throughput technologies in order to create comprehensive datasets [START_REF] Cox | Quantitative, high-resolution proteomics for data-driven systems biology[END_REF]. [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF] integrated proteome with transcriptome and metabolome data for the investigation of responses of E. coli cells to environmental and genetic perturbations.

Although proteomics can determine the absolute amount of each of the proteins in a sample or their relative change between several conditions, some major limitations remain. Efforts have to be made to increase its throughput compared to other large-scale technologies, by reducing measuring time on improved instruments [START_REF] Picotti | High-throughput generation of selected reaction-monitoring assays for proteins and proteomes[END_REF][START_REF] Reiter | mProphet: automated data processing and statistical validation for large-scale SRM experiments[END_REF]. Furthermore, special attention needs to be given to the development of simplified sample preparation protocols and computational analysis software.

As for high-throughput techniques, protein extraction from bacterial samples and the long processing times do not allow direct in vivo and real-time measurements [START_REF] Picotti | Selected reaction monitoring-based proteomics: workflows, potential, pitfalls and future directions[END_REF].

RT-qPCR

Real-time polymerase chain reaction (RT-PCR) is a recent technology developed for molecular biology and medicine, based on DNA labeling dyes to quantify changes in RNA expression levels. Quantitative RT-PCR has become one of the most popular methods for the analysis of gene expression and verification of microarray results.

While high-throughput microarray analysis allows large-scale analysis of gene expression profiles, the reverse transcription (RT) followed by the polymerase chain reaction (PCR) are often used to validate their findings [START_REF] Bustin | Pitfalls of quantitative real-time reverse-transcription polymerase chain reaction[END_REF][START_REF] Bustin | The need for transparency and good practices in the qPCR literature[END_REF]. In addition, relative quantification or absolute quantification of gene expression compared to standards that are run in parallel can be performed.

RT-qPCR protocols consist of several steps. First, RNA is isolated from sample cells, then mRNA is reverse-transcribed to cDNA. Second, the synthesized cDNA is amplified using specific PCR primers for the gene of interest. The PCR reaction also contains a fluorophore that specifically binds to double-stranded DNA, which allows

Experimental techniques for measuring gene expression

real-time monitoring of the amplification reaction. The fluorophores that are currently available for quantitative PCR can be broadly classified as non-specific and sequencespecific DNA-associating dyes [START_REF] Mackay | Real-Time PCR in microbiology: From diagnosis to characterization[END_REF][START_REF] Saunders | Real-Time PCR: Advanced Technologies and Applications[END_REF]. The most commonly used reporter is the SYBR Green dye, a non-specific probe. Once it is bound to the double-stranded DNA, it emits a fluorescent signal increasing with the accumulation of PCR products. The other specific reporters such as TaqMan rely on Fluorescence Resonance Energy Transfer (FRET) from the dye molecule to the quencher, both attached to a specific oligoprobe [START_REF] Holden | Quantitative Real-Time PCR: fluorescent probe options and issues[END_REF]. When the oligoprobe hybridizes with its template DNA, the fluorophores (dye molecule and quencher) are released and fluorescence emissions can be detected.

In conclusion, RT PCR has been mostly used for microarray data validation, but also for the absolute and relative quantification of the number of plasmid copies in [START_REF] Lee | Absolute and relative QPCR quantification of plasmid copy number in Escherichia coli[END_REF]. However, the method uses cell lysate and does not allow in vivo measurements.

Moreover, the technology is not easily parallelized, therefore it is only used for detailed study and validation of results. In addition, preparation of mRNA involves additional steps, may lead to the loss of some initial mRNA, and it is more difficult to assess the quality of the final product [START_REF] Saunders | Real-Time PCR: Advanced Technologies and Applications[END_REF].

Reporter genes

Current reporter gene technologies, based on Green Fluorescent Proteins (GFPs) [START_REF] Southward | The dynamic microbe: Green fluorescent protein brings bacteria to life[END_REF] and other fluorescent and luminescent reporter proteins, provide an excellent means to measure gene expression in vivo and in real time, in contrast to the other techniques presented in this chapter. The underlying principle of the technology is the fusion of the promoter region and possibly (part of) the coding region of a gene of interest and a reporter gene. The expression of the reporter gene generates a visible fluorescence or luminescence signal that is proportional to the actual number of fluorescent molecules [START_REF] Rosenfeld | A fluctuation method to quantify in vivo fluorescence data[END_REF], it is easy to detect and reflects the expression of a gene of interest when its promoter is activated. Reporter gene constructions have enabled the real-time tracking of gene expression dynamics in single cells, which provides valuable information about the functioning of cells [START_REF] Golding | Real-time kinetics of gene activity in individual bacteria[END_REF][START_REF] Longo | Dynamics of single-cell gene expression[END_REF]. Blue, yellow, cyan, and red fluorescent proteins have been engineered, allowing the expression of several genes to be studied in parallel, in the 2. STATE OF THE ART same cell. For example, three fluorescent reporter genes (encoding cyan, yellow and red fluorescent protein) were used to investigate how noise is transmitted in a gene network [START_REF] Pedraza | Noise propagation in gene networks[END_REF].

Reporter constructs can be located on plasmids or on the genome, depending on the

Reporter gene

Promoter region

Reporter protein

Figure 2.1: Reporter gene systems. The promoter region of a gene of interest is fused with a reporter gene. When the promoter is activated, the reporter gene is transcribed and produces the reporter protein as a direct quantification of the strength of the promoter.

The reporter system can be introduced in cells on a plasmid or directly integrated into the chromosome.

problem studied (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]. Plasmidic reporters have the advantage of being easy to construct and generating a strong signal, compared to chromosomal constructions, but their copy number may change with the experimental conditions (Lin-Chao and [START_REF] Lin-Chao | Effect of the bacterial growth rate on replication control of plasmid pBR322 in Escherichia coli[END_REF], thus introducing biases in the interpretation of the data. By measuring the plasmid copy number in a cell with qRT-PCR technology, one can correct this bias (Berthoumieux et al., 2013b). In addition, protein concentrations can be reconstructed from fluorescent reporter gene measurements by means of existing data processing methods (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF], when information about proteins half-lives is available.

Many studies use fluorescent reporter gene data to quantify gene expression in E. coli .

Inference of gene regulatory networks

For example, [START_REF] Ronen | Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics[END_REF] used fluorescent reporter genes to investigate its SOS DNA repair system. [START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF][START_REF] Kalir | A coherent feed-forward loop with a SUM input function prolongs flagella expression in Escherichia coli[END_REF] analyzed the dynamics of the gene network regulating flagella motif. To create a tool for accurate, high-resolution analysis of transcription networks, [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF] constructed a library of transcriptional fusions of gfp to the intergenic regions containing promoters, in E. coli K12 strain, on a low-copy plasmid. Dynamic measurements have been obtained using this library on a genomic scale, in a diauxic shift experiment [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF].

Inference of gene regulatory networks

Existing methods for gene network inference can be classified in different ways, depending on the criterion used for the classification. Several review papers explain the principles of the different inference methods. In [START_REF] Faith | Reverse-engineering transcription control networks[END_REF], inference methods are divided into 'physical interaction' approaches, that aim at identifying interactions among transcription factors and their target genes, and 'influence interaction' approaches, that try to relate the expression of a gene to the expression of other genes. [START_REF] Stelling | Mathematical models in microbial systems biology[END_REF] and [START_REF] Doyle Iii | Systems interface biology[END_REF] classify inference methods based on the type of models they pertain, distinguishing between interaction-based, constraint-based and mechanism-based models. A relevant classification of modelling and estimation approaches is provided in [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]. Discussions of modelling methods and details on the computational aspects can be also found in [START_REF] Beer | Predicting Gene Expression from Sequence[END_REF], [START_REF] Ambesi | Computational Biology and Drug Discovery: From Single-Target to Network Drugs[END_REF] and [START_REF] Markowetz | Inferring cellular networks-a review[END_REF]. A mixed classification of algorithms based on methodological approach, modelling context and performance on different inference problems results from the DREAM challenge (Dialogue for Reverse Engineering Assessments and Methods), a large and cooperative effort toward the assessment of inference performance [START_REF] Greenfield | DREAM4: Combining genetic and dynamic information to identify biological networks and dynamical models[END_REF][START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF][START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF] in the form of a competitive game among participant methods.

In what follows we will review gene network inference algorithms associated with different interaction modelling formalisms, mostly following the recent review by [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]. In every case we will discuss the main approaches proposed in the literature along with their advantages and limitations.

STATE OF THE ART

Inference of gene clusters

A first approach to gene network inference is to try group genes that may be functionally related. Genes whose expression appears to be altered in a coordinated manner in response to one or more experimental perturbations are grouped together and considered to be related to the cell function that has been probed experimentally. Typically in this context, gene expression screening is performed by way of microarray experiments, and functional gene grouping is operated by way of statistical methods called clustering.

Clustering relies on the concept of similarity among expression patterns or profiles [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]. As a similarity metrics, a correlation coefficient, most commonly the Pearson coefficient, is used:

r ij = M k=1 (x i (k)x j (k)) M k=1 (x 2 i (k)) M k=1 (x 2 j (k)) (2.1)
where x i and x j are gene expression measurements taken in M different conditions and r ij is the pairwise correlation coefficient computed between gene i and gene j. All pairwise correlation coefficients r ij for all possible gene pairs ij are computed for a set of n profiles. If the expression patterns of two genes are perfectly correlated (i.e., they are identical up to shifting and scaling) then r ij = 1; in the opposite situation, when the variables are linearly independent, r ij = 0. Correspondingly, gene pairs with large enough correlation coefficients are deemed to be functionally related (e.g.

above a suitable threshold between 0 and 1). Based on this principle, suitable methods for grouping genes into similarity clusters have been developed [START_REF] Amato | A multi-step approach to time series analysis and gene expression clustering[END_REF][START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]. Gene clusters obtained in this way can be seen as fully connected subgraphs of the graph with all genes as the graph nodes. Links among nodes within a cluster are thus functional relationships and, although no direct mutual regulation among clustered genes is implied by this grouping, clustered genes are often presumed to undergo some form of mutual interaction. A common subsequent analysis step is to annotate each cluster with a functional category representative of that cluster and use this categorization for further inference [START_REF] Guthke | Dynamic network reconstruction from gene expression data applied to immune response during bacterial infection[END_REF].

Inference of gene regulatory networks

A common issue in gene clustering is the choice of the policy for associating genes to one cluster or another. In particular, the number of clusters is a priori unknown and can be chosen automatically based on data or manually, depending on the clustering algorithm used [START_REF] Amato | A multi-step approach to time series analysis and gene expression clustering[END_REF][START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF].

Similarity measures alternative to (2.1) have also been proposed, and the choice has important effects on the clustering results. However, validation of results still remains an opened issue [START_REF] Allison | Microarray data analysis: from disarray to consolidation and consensus[END_REF][START_REF] Handl | Computational cluster validation in post-genomic data analysis[END_REF]. Most current clustering algorithms do not provide estimates of the significance of the results returned. The validation of clustering results is therefore often based on a manual and subjective exploration process, such as visual inspection and prior biological knowledge to select what is considered the most "appropriate" result. Recently, clustering has been applied

to metabolic pathways analysis [START_REF] Milone | Improving clustering with metabolic pathway data[END_REF], incorporating prior knowledge into the cluster formation itself and show important improvements in the convergence and performance.

Most clustering methods developed for gene networks [START_REF] Liang | Genetic network inference: from co-expression clustering to reverse engineering[END_REF][START_REF] Stuart | A gene-coexpression network for global discovery of conserved genetic modules[END_REF] are derived from hierarchical clustering [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]. (2007)) some known interactions were indeed recovered. Newer, sophisticated methods such as distance correlation have theoretical advantages over Pearson coefficient [START_REF] Székely | Brownian distance covariance[END_REF] and have been tested on protein networks [START_REF] Roy | Detection of long-range concerted motions in protein by a distance covariance[END_REF].

Recently, methods for clustering have been widely coupled to more complex inference algorithms to reduce the dimensionality of the search space before network inference (cMonkey, [START_REF] Reiss | Integrated biclustering of heterogeneous genome-wide datasets for the inference of global regulatory networks[END_REF]).

Inference of interaction graphs

Clustering approaches based on the Pearson coefficient may be effective for linear correlations, but their performance decreases for nonlinear systems (Villaverde and Banga,

STATE OF THE ART

X 2 X 3 X 1 X 4 X 5 X 2 X 3 X 1 X 4 X 5 X 2 X 3 X 1 X 4 X 5 AND X 2 X 3 X 1 X 4 X 5 X 2 X 3 X 1 X 4 X 5
Time series

Gene cluster

Information-theoretic approach

Boolean network

Bayesian network

Ordinary differential equations

Gene expression

Steady-state

} , min{ 0 , 0 34 13 14 12 15 MI MI MI MI MI    ) ( ) ( ) 1 ( 3 2 1 t x t x t x    ) , | ( ) , , | ( 3 2 1 4 3 2 1 x x x P x x x x P  ) , , , (
) ( 3 2 1 1 1  x x x g dt t dx      2 2 2 1 2 1 12 x x x x r

INFERENCE

Figure 2.2: Approaches for inferring gene regulatory networks. In gene clusters, all dependencies between genes are typically assessed by Pearson correlations. In information-theoretic networks, MI is 0 for statistically independent variables and DPI is used to select direct regulatory interactions. In Boolean networks, the state of a gene is computed as a simple Boolean rule from the activities of other genes. Bayesian networks employ probability distributions to determine regulatory effects between genes. In the case of ODE models, the activity of one gene is computed as a function (g) of the level of its regulators. (Inspired by [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF] 2014). In addition, being the focus on functional relationships, no distinction is usually made between direct and indirect interactions. A related problem is that of inferring graphs of interactions, where nodes are genes but edges represent actual regulatory

Inference of gene regulatory networks

effects between gene pairs. To discriminate between direct and indirect interactions, rather than looking at Pearson correlation (Eq. 2.1), information-theoretic concepts such as Mutual Information (MI, [START_REF] Shannon | A mathematical theory of communication[END_REF]) can be employed instead. Based on this, in general agreement with literature terminology, we will discuss inference of interaction graphs in terms of methods based on an information-theoretic approach.

The underlying principle of information-theoretic approaches is the concept of entropy, describing the uncertainty of a single random variable:

H i = - x i ∈χ p(x i )log(p(x i )) (2.2)
which captures the a priori variability of the expression of the gene, and more generally, the joint entropy

H ij = - x i ∈χ x j ∈χ p(x i , x j )log(p(x i , x j )) (2.3)
which quantifies the variability of the random variables involved. Mutual information, M I ij , between gene i and gene j is computed as:

M I ij = H i + H j -H ij (2.4)
It is a measure of dependencies between variables, the higher the value, the stronger the mutual dependency. If two variables are statistically independent, their joint en-

tropy is H ij = H i + H j (i.e. p(x i x j ) = p(x i )p(x j
)) and the mutual information is zero.

Conversely, full dependency corresponds to H ij = 0. In practice, these quantities can be computed e.g. from microarray gene expression profile data [START_REF] Butte | Mutual information relevance networks: functional genomic clustering using pairwise entropy measurements[END_REF], see further below. Similar to the correlation approach, choosing a MI threshold provides a mechanism to identify potential regulatory interactions, the smaller the threshold, the higher the hit rate at the price of more probable false hits.

Different from the correlation approach, use of the so-called data processing inequality allows one to prune many indirect interactions among genes. The data processing inequality states that if genes (i, k) interact indirectly through j, and no alternative path exists between genes (i, k), then M I ik ≤ min(M I ij , M I jk ). Thus, when the data processing inequality is verified, indirect interactions can be excluded from the set of estimated interactions in favor of direct interactions. A well known example of inference method based on the DPI is ARACNE [START_REF] Basso | Reverse engineering of regulatory networks in human B cells[END_REF][START_REF] Margolin | ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF].

In ARACNE, M I ij is estimated from data for all pairs of observed genes i and j by

STATE OF THE ART

the use of a Gaussian kernel [START_REF] Steuer | The mutual information: detecting and evaluating dependencies between variables[END_REF]. Based on Montecarlo simulations, a threshold M I 0 corresponding to a suitable p-value for testing the independence of Gaussian variables is determined. All gene pairs with mutual information below this threshold are considered independent. Putative interactions are then checked for false positives, i.e. dependence relationships that are not direct are sought based on the DPI.

This allows elimination of indirect interactions but may also eliminate direct interactions [START_REF] Margolin | ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF], depending on the data and the tuning of certain algorithm parameters. Automatic choice of MI and DPI threshold parameters suggested along with ARACNE gives a good sensitivity-accuracy compromise [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF].

In order to further improve the accuracy vs. sensitivity tradeoff, other approaches have been investigated, leading e.g. to the design of minimum redundancy [START_REF] Meyer | Information-theoretic inference of large transcriptional regulatory networks[END_REF], entropy reduction [START_REF] Samoilov | Reconstruction and functional analysis of general chemical reactions and reaction networks[END_REF][START_REF] Villaverde | Reverse Engineering Cellular Networks with Information Theoretic Methods[END_REF] and continuous threeway mutual information [START_REF] Luo | Learning transcriptional regulatory networks from high throughput gene expression data using continuous three-way mutual information[END_REF] methods.

Since M is symmetric (M ij = M ji ) and mutual information does not provide any information about directionality of regulation, the interaction network is still reconstructed in the form of an undirected graph G, where edges represent statistical dependencies observed in the data, but do not carry information about causality of regulation. In order to establish causal relationships from the inferred associations between interacting nodes, i.e. discriminate between regulatory and target genes, additional information is necessary. Toward this aim, one example is Context Likelihood of Relatedness (CLR) algorithm [START_REF] Faith | Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF] where a distinction is made between the roles of transcription factors and target genes. This algorithm includes an adaptive background correction to reduce false corrections and indirect influences. CLR calculates the statistical likelihood of each MI value and then compares the MI of a transcription factor/gene pair to the background distribution of all possible transcription factor/gene pairs that include either the transcription factor or its target. Retaining only regulations whose associated MI is significantly higher than the background values, apparent interactions stemming from e.g. activity of one transcription factor weakly varying with the expression of several genes, or expression of one gene weakly varying with the activity of several transcription factors, are filtered out of the reconstructed network.

In practice, the ability to estimate the MI between genes from experimental data, and hence the practicality of the methods discussed above, depends on the quality of the data, in particular, on the relevance of certain statistical independence assumptions.

Inference of gene regulatory networks

Estimation of MI can be performed from multiple independent steady-state datasets as well as from dynamic data, as long as the sampling time is long enough to consider subsequent measurements statistically independent [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]. Computation of MI indices then involves the probability mass function p(x), which is typically unknown in reverse-engineering problems and needs to be reconstructed from the data. This is often done by histogramming, i.e., partitioning the data into equally sized bins and then counting the frequency of appearence of the data in every bin [START_REF] Steuer | The mutual information: detecting and evaluating dependencies between variables[END_REF]. This step is important and may affect the gain in performance from

Pearson correlation-based methods, quantifying only linear dependencies, and mutual information-based methods [START_REF] Steuer | The mutual information: detecting and evaluating dependencies between variables[END_REF].

In real-world applications, ARACNE has been shown to perform well on large and medium size steady-state datasets (human Bcells and S. cerevisiae, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]), yielding results similar to clustering algorithms. On small time-varying datasets, ARACNE had a poor performance, as expected because of its requirement of statistically independent time-points [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]. Therefore, a new version of ARACNE better suited for time-course data has been recently developed [START_REF] Zoppoli | TimeDelay-ARACNE: Reverse engineering of gene networks from time-course data by an information theoretic approach[END_REF].

Using a compendium of microarray expression profiles in E. coli, CLR was used in [START_REF] Faith | Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF] not only to reconfirm known regulations, but also to discover several novel interactions, some of which were validated experimentally. Still, direction of regulation is well defined only for interactions discovered between one transcription factor and one target gene, while it is undefined for interactions found between transcription factors.

In combination with other methods, CLR was reported to be an effective inference algorithm for the DREAM4 100-gene in silico network inference challenge [START_REF] Greenfield | DREAM4: Combining genetic and dynamic information to identify biological networks and dynamical models[END_REF].

Thus, despite the variously successful applications and the development of methods such as CLR, the main challenge of the approaches developed in this framework remains that of establishing the causality of interactions. In the next section, a framework that allows for a natural treatment of the direction of regulations is discussed.

STATE OF THE ART

Inference of Boolean networks

Boolean networks are discrete network models (or logical networks). Boolean networks were first described by [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF] and were later on recognized as a natural framework for gene regulation modelling [START_REF] Bornholdt | Boolean network models of cellular regulation: prospects and limitations[END_REF][START_REF] Kauffman | Random Boolean network models and the yeast transcriptional network[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF]. The state of a gene (nodes in the network) can be approximated or described by a Boolean variable x i ∈ {0, 1}. The gene can be "inactive"(0) or "active"(1) (de [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]. Boolean networks are usually represented as directed graphs, where the edges (interactions between nodes) are represented by activation/inhibition Boolean functions. Through these functions, the state of a gene is determined on the basis of the states of its parent genes by applying basic Boolean operations (AND, OR, NOT).

Boolean networks can capture the dynamics of a regulatory system on a discrete time grid as follows:

x i (t + 1) = g i (x 1 (t), ...x k (t)) (2.5)
where x i (t), the state of the gene i at time t changes to a state, x i (t + 1) at time t + 1 following Boolean function g i . Each gene can be active or inactive, so the state space consists of 2 N states. The boolean function uses states for k nodes (regulators), so the number of possible functions is 2 2 k (de [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]. Although rapidly increasing with the number of genes in the network, the state space and the number of interaction functions is finite. Therefore, Boolean networks evolve towards a steady state or a cycle of states, called attractor [START_REF] Klipp | Growth-rate-dependent partitioning of RNA polymerases in bacteria[END_REF] et al., 2000;[START_REF] Shmulevitch | Probabilistic Boolean Networks: a rule-based uncertainty model for gene regulatory networks[END_REF] have been developed as an alternative to classic deterministic Boolean networks to deal with data uncertainty, allowing more than one possible state transition Boolean function to describe gene interaction from noisy expression data. Nevertheless, Boolean networks provide a good qualitative interpretation of gene regulation and can be used to simulate gene regulatory functions [START_REF] Hecker | Gene regulatory network inference: data integration in dynamic models-a review[END_REF].

Inference of Bayesian networks

Bayesian Networks (BN) are probabilistic graph-based methods characterizing the expression of every gene i in a regulatory network by a random variable X i . The interactions among these are represented as joint probability distributions P (X 1 , ..., X N ) and encoded in the structure of a directed acyclic graph G, whose nodes are the random variables X i . The joint probability density is expressed as a product of conditional probabilities by applying Bayes' theorem:

P (A, B) = P (B A) * P (A) = P (A B) * P (B).
This allows one to write

P (X 1 , ..., X N ) = N i=1 P (X i = x i X j 1 = x j 1 , ..., X jp = x jp ).
(2.6)

The p genes (with p generally depending on i) that appear in the conditional probability for gene i represent the actual regulators of gene i, i.e. the directed edges of the associated graph G (also called the parents of node i). A crucial assumption is that the model obeys a Markov property stating that each variable X i is conditionally independent of its non-descendants given its parents. If this causality condition does not hold, similar to information-theoretic approaches, the edges of the BN graph no longer represent direct causal interactions but only statistical dependencies [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF].

Inferring a Bayesian network model means finding the directed acyclic graph G that best explains a gene expression dataset, in the sense of maximizing a probabilistic scoring function related to Eq. (2.6) over the candidate network topologies [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF][START_REF] Nir | Inferring cellular networks using probabilistic graphical models[END_REF][START_REF] Villaverde | Reverse engineering and identification in systems biology: strategies, perspectives and challenges[END_REF]. Existing methods typically explore the multidimensional space of possible graphs (i.e., regulatory networks) exhaustively [START_REF] Faith | Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF][START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF].

STATE OF THE ART

In reverse engineering applications, neither prior nor posterior probabilities are known [START_REF] Villaverde | Reverse engineering and identification in systems biology: strategies, perspectives and challenges[END_REF]. Markov Chain Monte Carlo (MCMC) techniques [START_REF] Gelfand | Sampling-Based Approaches to Calculating Marginal Densities[END_REF][START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF] are often used to evaluate Eq. (2.6) and other relevant probabilities. Still, the search for the best BN model is an NP-hard problem [START_REF] Chickering | Learning from data: artificial intelligence and statistics[END_REF]. Hence, heuristic solutions are often implemented. In the context of the DREAM project, MCMC BN inference was found to be computationally very costly, making this approach better fitted for smaller networks [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF].

In addition, the BN learning problem is usually underdetermined. Several strategies exist to deal with this underdetermination, such as model averaging, bootstrapping or if available, adding a priori knowledge to select the most likely model structure [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF][START_REF] Villaverde | Reverse engineering and identification in systems biology: strategies, perspectives and challenges[END_REF]. BN inference is well suited to integrate heterogeneous datasets [START_REF] Klipp | Systems biology in practice: concepts, implementation and application[END_REF], which makes BN modelling appealing when compared to the network inference methods reviewed above.

However, the main limitation of Bayesian network models is that they cannot contain cycles and thus they cannot represent cellular network feedback loops [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF][START_REF] Doyle | Systems interface biology[END_REF]. To overcome this limitation, Dynamic Bayesian networks were developed along with inference methods using time-series data [START_REF] Yu | Advances to Bayesian network inference for generating causal networks from observational biological data[END_REF].

A well-known BN inference algorithm is Banjo [START_REF] Yu | Advances to Bayesian network inference for generating causal networks from observational biological data[END_REF]. This algorithm addresses both static and dynamic BNs and hence can infer gene networks from steadystate and dynamic data, based on heuristic approaches for the exploration of the possible network topologies. The algorithm was tested on simulated and real microarray data [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]. Banjo showed poor performance when applied on steadystate data from a limited number of experiments (E. coli steady-state dataset, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]) when compared to interaction graph inference via ARACNE. In presence of more data from a larger number of steady-state experiments (S. cerevisiae, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]), inference accuracy of Banjo proved to be considerably better (S. cerevisiae, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]) although the number of correct inferred interactions remained limited. Unfortuntately, applicability of Banjo is limited for increasing size of the dataset (number of experiments) due to computational complexity (HumanBcells and S. cerevisiae, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]). In the case of dynamic data analysis, relative to the number of network genes, analysis based on both synthetic 2.2 Inference of gene regulatory networks and experimental data (E. coli steady-state, see [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]) showed that Banjo [START_REF] Yu | Advances to Bayesian network inference for generating causal networks from observational biological data[END_REF] requires a large number of experiments for the estimation of the necessary probability distributions.

Inference of Ordinary Differential Equation models

We have seen so far inference of models describing qualitative relationships among genes, whether in the form of statistical relationships or in the form of regulatory logics. What all these frameworks are missing is quantitative information concerning the strength with which genes regulate each other, and how this reflects into quantitative time-course dynamics. In this section we look at Ordinary Differential Equation (ODE)-type modelling, which will also constitute the modelling framework of choice in the following chapters.

In the context of gene networks, ODE models describe the dynamics underlying timecourse gene expression in terms of a time-varying network state. The state generally represents concentrations of gene products (proteins, and sometimes corresponding mRNAs) and transcription factors that control the regulation of the genes of interest. Thus, this approach provides a more detailed and complex representation of the functioning of the biological network [START_REF] Villaverde | Reverse Engineering Cellular Networks with Information Theoretic Methods[END_REF] and is better suited for the analysis of and inference from time-course quantitative (population-average) gene expression data.

In their simplest form, ODE models of gene regulation have the representation

dx i (t) dt = g i (x 1 , ..., x N , θ i , u(t)) -λ i x i (t), i = 1, . . . , N (2.7)
where x i (t) is the time-varying concentration of the product of gene i, λ i is a degradation rate, g i is the synthesis rate function, θ i is a vector of characteristic parameters and u represents a network input or perturbation (chemical treatments, genetic modifications, etc.) for a network with N genes [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF].

From a gene interaction viewpoint, these models encode the regulatory effects of every gene on gene i into the analytical form of g i , thus enabling a link with the previously discussed model classes. While the regulatory structure can in principle be captured by suitably defined parameters θ i , it is generally more appropriate to keep the concepts of structure and continuous-type parameters distinct. The initial state of the system, x i (t 0 ) with i = 1, . . . , N , often is a steady-state solution of Eq. (2.7) (for an

STATE OF THE ART

appropriate constant pre-experimental input) corresponding to the system being at a dynamical equilibrium at the beginning of the experiment (de [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]. Inference of ODE models most often refers to the problem of estimating parameters θ i from timecourse or multiple steady-state experimental measurements of concentrations x i from one or more system perturbations u, for a given structure of functions g i . However, much work has also been dedicated to reconstructing or selecting the structure of the g i , whence notably the set of regulatory interactions among the network genes, from within suitable model families (Porreca et al., 2010a).

Unlike other approaches discussed in this chapter, in accordance with its deterministic nature, ODE modelling usually does not rely on a statistical characterization of gene regulatory interactions.

Model inference requires first of all the choice of an appropriate mathematical representation of g i , whether a model family or a fixed model structure, usually from the class of low-order polynomials or a combination of Hill functions (see [START_REF] Aracena | Maximum number of fixed points in regulatory Boolean networks[END_REF]; Given a model structure, the identification of nonlinear ODE parameters from gene expression data is a challenging problem, since the ODE system often does not have an explicit solution (de [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]. In general, without suitable constraints, there are multiple solutions, i.e. the ODE system is not uniquely identifiable from data. A common distinction is made between structural identifiability and practical identifiability.

Despite lack of agreement in the literature, the first generally refers to impossibility to distinguish different parameter values for the given model and observed outputs, no matter the abundance of the data, while the second refers to inherent limitations in estimation accuracy due to the quality of the data (Berthoumieux et al., 2013a).

Indeed, time-course gene expression data are usually sparse and associated with large noise [START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF], which makes estimation of structurally identifiable parameters quite uncertain.

Parameter estimation is generally expressed as an optimization problem, where the objective functions minimized quantify the distance between the observations and the

Inference of gene regulatory networks

model-predicted values. Thus, model inference is performed using optimization techniques [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]. Convex optimization, when applicable, provides a unique minimum and scales well with the dimension of the problem (number of unknown parameters, etc.) [START_REF] Boyd | Convex optimization[END_REF] 2012)). However, their computational complexity increases rapidly with the problem size [START_REF] Miró | Deterministic global optimization algorithm based on outer approximation for the parameter estimation of nonlinear dynamic biological systems[END_REF].

Models inferred from data carry information of the interactions among genes and hence provide implicitly a signed graph of regulatory interactions. Different from previously discussed methods for network inference, however, ODE models (Eq. 2.7) with parameters inferred from data can then be used to predict quantitatively the time response of the network to different internal perturbations (e.g. gene knock-out or over-expression) and external stimuli (environmental changes).

Among The network is modeled as a system of linear ODEs (de [START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF] expressing the synthesis rate of every transcript as a linear function of the concentration x i of all other cell transcripts and a network perturbation u. For one experiment with given perturbation u, at measurement time t k , with k = 1, . . . , M , the system satisfies

ẋi (t k ) = N j=1 a ij x j (t k ) + b i u(t k ), (2.8)
where, for N genes, i = 1, ..., N . The coefficients a ij and b i (collectively captured by θ i in Equation 2.7) quantify the effect of gene j and perturbation u, respectively, on gene i. If u is constant and the system has reached steady-state, then ẋi (t k ) = 0 and, for all 2. STATE OF THE ART i, Equation 2.8 simplifies to

N j=1 a ij x j (t k ) = -b i u(t k ) (2.9)
From steady-state gene expression data (M = 1 and measurements of x i for different constant values of u), the NIR algorithm [START_REF] Gardner | Inferring genetic networks and identifying compound mode of action via expression profiling[END_REF] computes the edges a ij by solving the multiple linear regression (2.9) assuming a priori information about genes that have been directly perturbed (i.e. terms b i u known). The user can choose the maximum number of regulators per gene (i.e., the number of edges to a node).

After solving linear regression (2.9), the algorithm returns the estimated matrix of interaction strenghts a ij . If the noise in the data is small, this method does not require large datasets. It has been tested in [START_REF] Bansal | Inference of gene regulatory networks and compound mode of action from time course gene expression profiles[END_REF] showing good performance compared to ARACNE or Banjo even when a reduced number of experiments are available.

Similarly, MNI algorithm (di [START_REF] Bernardo | Chemogenomic profiling on a genome-wide scale using reverseengineered gene networks[END_REF] is based on Equation 2.9 and uses steady-state data. However, each microarray experiment can result from any kind of perturbation and knowledge about b i u is not necessary. First, MNI computes the a ij from the gene expression data D and determines a model of the regulatory interactions between genes. Then a test dataset {x d 1 , ..., x d N } representing the perturbed expression of the genes is used to compute b i u from Equation 2.9, with u a simple constant. The network model initially identified in the algorithm (trained on the dataset D) is used as a filter to predict a better model from the test perturbation data. A b i = 0 quantifies that gene i is directly affected by the perturbation. The algorithm returns a ranked list of genes, where most likely targets of the perturbation have high values of b i .

The TSNI algorithm [START_REF] Bansal | Inference of gene regulatory networks and compound mode of action from time course gene expression profiles[END_REF] relies on dynamical time-series data and identifies both the network structure a ij and the targets of perturbation, that is the b i , by solving in this case the linear regression corresponding to a discrete-time version of Equation 2.8. The algorithm assumes that a single perturbation experiment is performed and M time points following the perturbation are measured, in contrast to M different steady-state conditions considered for NIR and MNI (it can, however, be easily generalized to multiple dynamical experiments). The algorithm is capable of correctly inferring the structure (a ij ) and the targets of perturbation (b i ) of small gene networks. In larger networks instead, its performance in recovering the structure is not 2.2 Inference of gene regulatory networks very good, probably because the network is not fully observable and one perturbation experiment does not yield sufficient information for the inference (SOS E. coli network, see [START_REF] Bansal | Inference of gene regulatory networks and compound mode of action from time course gene expression profiles[END_REF][START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]).

Many other ODE-based algorithms have been proposed in the literature. Inferelator [START_REF] Bonneau | The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo[END_REF]) is an inference algorithm for deriving genome-wide transcriptional regulatory interactions, and has been applied to predict a large portion of the regulatory network of the archaeon Halobacterium NRC-1. The algorithm infers regulatory interactions for genes and gene clusters from mRNA or protein expression levels and uses standard regression and l 1 shrinkage techniques to select models for the expression of a gene or cluster of genes as a function of the levels of their regulators. In [START_REF] Bonneau | The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo[END_REF], many novel gene interactions were predicted, and in several cases the inferred regulatory interactions were validated by experimental tests. The Inferelator was also able to predict mRNA levels of 80% of the genes in the genome over new experimental conditions in Halobacterium salinarium [START_REF] Bonneau | A predictive model for transcriptional control of physiology in a free living cell[END_REF]. [START_REF] Greenfield | DREAM4: Combining genetic and dynamic information to identify biological networks and dynamical models[END_REF] demonstrate complementarity between this method and the mutual information CLR (Section 2.2.2) algorithm. Based on application to in silico time-series data of the DREAM4 competition, their combined use significantly improves the ability of selecting valid regulatory interactions compared to both methods alone.

Moreover, the duo is able to accurately predict the response of the system to new conditions (new double knock-out perturbations).

Inference of Boolean-like models and the sign pattern analysis method

In the spirit of mixed ODE-boolean modeling, an original method of network inference has been proposed in Porreca et al. (2010a), explicitly based on the idea of transcription rate functions encoding the logics of regulation of the target genes. The resulting ODE models are referred to as "Boolean-like" models. The proposed inference algorithm tackles the identification of both structure and parameters of kinetic models of gene regulatory networks from time-course gene expression data. A modeling framework is considered where the dynamic equations are described in terms of a class of gene activation rules known as unate functions [START_REF] Aracena | Maximum number of fixed points in regulatory Boolean networks[END_REF][START_REF] Comet | On circuit functionality in boolean networks[END_REF].

These functions reflect interactions where each gene is exclusively either an activator or an inhibitor for the expression of any given target gene (though a regulator may 2. STATE OF THE ART be an activator and a repressor of distinct target genes). According to [START_REF] Grefenstette | An analysis of the class of gene regulatory functions implied by a biochemical model[END_REF], unate functions provide biologically realistic dynamics for gene networks.

The majority of the known gene activation rules [START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF][START_REF] Nikolajewa | Boolean networks with biologically relevant rules show ordered behavior[END_REF] are modeled by nested canalizing functions [START_REF] Jarrah | Nested Canalyzing, Unate Cascade, and Polynomial Functions[END_REF], a class of unate functions. Unate function modelling and analysis of biochemical networks has been discussed also in [START_REF] Murrugarra | Regulatory patterns in molecular interaction networks[END_REF][START_REF] Raeymaekers | Dynamics of Boolean networks controlled by biologically meaningful functions[END_REF].

In Porreca et al. (2010a), the properties of these functions are exploited in order to develop a two-step identification algorithm. Boolean-like ODE models are used to describe the evolution of the product of gene i (x i ):

ẋi = g i (x) -γ i (x)
(2.10) where x = (x 1 , ..., x N ), N is the number of genes in the network, while g i (x) ≥ 0 and γ i (x) ≥ 0 are the synthesis and the degradation rates of the product of gene i.

The nonlinear model for the synthesis rate is

g i (x) = k 0,i + k 1,i b i (x) (2.11)
where k 0,i ∈ R + and k 1,i ∈ R + are constants and b i (x) : R n + [0, 1] quantifies the regulatory effects of the gene products in the network on the expression of gene i by algebraic combinations of Hill activation or repression functions [START_REF] Keller | Model genetic circuits encoding autoregulatory transcription factors[END_REF][START_REF] Yang | An analytical rate expression for the kinetics of gene transcription mediated by dimeric transcription factors[END_REF],

namely σ + (x j ) = x d j x d j + θ d , σ -(x j ) = 1 -σ + (x j ).
(2.12)

Due to the assumption of unate structure, every function b i , and hence every corresponding g i , is monotonically increasing or decreasing in every state variable x that is an effective regulator of target i, while it is independent of x j if the product of gene j does not regulate expression of gene i. These monotonicity properties can be captured by a sign pattern, i.e. an N -tuple p = (p 1 , ..., p N ) ∈ {-1, 0, 1} N (depending on i) where, for j = 1, ..., N p j is -1 if gene j acts as an inhibitor for gene i expression, it is 1 if gene j acts as an activator for gene i expression, and 0, if gene j has no effect on the expression of gene i.

The problem tackled by the two-stage algorithm is the reconstruction of every g i . To

this purpose, and unlike other algorithms in this section, here the authors assume measurements of time-varying protein concentrations and promoter activities. As measurements of x and corresponding g i (x) are considered to be available, the decay rate 2.2 Inference of gene regulatory networks γ i (x) can be ignored in the reconstruction of g i . In practice, when only concentration measurements or synthesis rate measurements are available, the dataset required by the algorithm can be completed from the available data provided knowledge of γ i , see Porreca et al. (2010b).

In the first step, the algorithm isolates families of consistent model structures, i.e.

families of so-called consistent sign patterns, by testing hypothetical sign patterns and rejecting those corresponding to monotonicity properties of g i that are falsified by the experimental data. This reduces considerably the number of plausible interactions. In addition, the family of consistent patterns can be arranged in a hierarchical fashion, and is fully characterized by a small set of minimal possible topologies of the network.

In the second step, quantitative identification of the networks returned by the first step is performed. By solving a nonlinear regression problem (estimation of θ, d, k 0,i and k 1,i and selection of the best specific model structure among those with a given sign pattern), models of minimal complexity explaining the data with sufficient accuracy are returned.

The method has been tested on an in silico network and on real data from synthetic network (IRMA, see [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverseengineering and modeling approaches[END_REF]) and compared to TSNI. The signed directed graphs inferred from IRMA by TSNI were less accurate than those of Porreca et al. (2010a), where an analysis of sensitivity to noise of network reconstruction performance was also provided based on in silico data.

Performance of existing inference algorithms, data requirements and perspective for novel developments are largely discussed by DREAM reports and in [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF]; [START_REF] Hecker | Gene regulatory network inference: data integration in dynamic models-a review[END_REF]; [START_REF] Villaverde | Reverse Engineering Cellular Networks with Information Theoretic Methods[END_REF]. Reliable inference from gene expression data still remains an open subject. One of the main conclusions is that, indeed, the performance of current network-inference methods is strongly dependent on the properties of the network that is being inferred and cannot be analyzed in isolation of the data that made it necessary. Employing modeling formalisms and algorithms that train on different features of the data and merging results seems to be a good way to improve the performance of inference. Another important point is that models and inference methods should be interpretable in terms of biological relevance of the results. For instance, while steady-state methods can assume protein and mRNA concentrations to be correlated, this can lead to spurious inference results from time-series datasets. However, the access to dynamic measurements for these chemical species is

STATE OF THE ART

not always easy and computational methods need to be developed in order to account for these issues. Next section addresses this point in a way that will be used for the results of this thesis.

Reporter gene data analysis 2.3.1 Measurement models

The reconstruction of biologically-relevant quantities from reporter gene data requires measurement models making explicit the relations between the concepts and the assumptions under which these relations hold (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]. Measurement models can describe the expression of the gene of interest in two steps (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]:

d dt m(t) = g(t) -(µ(t) + γ m ) m(t), m(0) = m 0 , (2.13)
d dt p(t) = κ p m(t) -(µ(t) + γ p ) p(t), p(0) = p 0 , (2.14)
where m(t), p(t) are the mRNA and protein concentrations, respectively, µ(t) is the time-varying growth rate, κ p is the protein synthesis rate constant, and γ m , γ p are the degradation constants of mRNA and protein, respectively. A similar measurement model can be written for the reporter protein:

d dt n(t) = g(t) -(µ(t) + γ n ) n(t), n(0) = n 0 , (2.15) d dt r(t) = κ r n(t) -(µ(t) + γ r ) r(t), r(0) = r 0 , (2.16)
with analogous meanings for the variables and parameters. By construction of the transcriptional fusions, the mRNA synthesis rates or promoter activities of the gene of interest and the reporter gene are equal. This promoter activity is denoted by g(t).

Two common assumptions make it possible to simplify the above models. First of all, typical mRNA half-lives in bacteria are on the order of a few minutes [START_REF] Bernstein | Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays[END_REF], whereas typical cell doubling times range from tens of minutes to hours [START_REF] Larrabee | The relative rates of protein synthesis and degradation in a growing culture of Escherichia coli[END_REF][START_REF] Mosteller | Metabolism of individual proteins in exponentially growing Escherichia coli[END_REF]. This motivates γ m , γ n µ(t). Second, the mRNA concentrations evolve on a much faster time-scale than the protein concentrations, so that the former can be assumed to be in quasi-steady state: dm(t)/dt = dn(t)/dt = 0.

Reporter gene data analysis

As a consequence, m(t) = g(t)/γ m and n(t) = g(t)/γ n , and the models of Eqs. 2.13-2.16 simplify to the following reduced models:

d dt p(t) = kp g(t) -(µ(t) + γ p ) p(t), p(0) = p 0 , (2.17) d dt r(t) = kr g(t) -(µ(t) + γ r ) r(t), r(0) = r 0 , (2.18)
with kp = κ p /γ m and kr = κ r /γ n . We define the synthesis rate of the reporter protein

f (t) = kr g(t). (2.19)
This quantity is proportional to the synthesis rate of the protein of interest, with

proportionality constant α = (κ r /κ p )(γ m /γ n ), i.e., f (t) = α kp g(t).
(2.20)

Therefore, if κ p = κ r (true for translational fusions) and γ m = γ n , then f (t) also equals the synthesis rate of the protein of interest. As explained in [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF], f (t) can be directly computed from the absorbance and fluorescence signals. The quantity is usually called promoter activity in the literature or more generally the activity of the gene, motivated by the fact that it is proportional to g(t). Promoter activity is also proportional to the mRNA concentration of the gene of interest. This simply follows from the fact that f (t) is proportional to kp g(t) and the latter expression equals k p m(t) by Eq. 2.14.

Constitutive promoters

One of the limitations of the above measurement model is that it assumes that k p , k r (the protein and reporter protein synthesis rate constants) are constants and do not depend on the time-varying activity of the ribosomes. The model also does not distinguish between the contributions of specific transcription regulators and the activity of RNA polymerase to the promoter activity g(t). In order to address these limitations, the measurement models can be easily generalized (Berthoumieux et al., 2013b)

proportionality constant (κ r /κ p )(γ m /γ n )).
If we consider a reporter gene with a constitutive promoter that has the same ribosomebinding site as the reporter of the gene of interest, following Eq. 2.22, we have:

g const (t) = k const m g global (t), (2.23)
and, correspondingly,

f const (t) = k const m kr (t) g global (t), (2.24)
Therefore, when measuring both f (t) (by means of the reporter of the gene of interest) and f const (t) (by means of the reporter of a constitutively expressed gene), global physiological effects due to the activity of the gene expression machinery and specific effects due to transcription factors and other regulators can be separated.

Data processing

As described in Section Experimental data of this chapter, in vivo and real-time gene expression profiles can be obtained by means of fluorescent reporter gene systems monitored in an automated, thermostated reader. The absorbance or the optical density measured at 600 nm quantifies the biomass, while the fluorescence signal emitted at 520 nm, when excited at 485 nm, is proportional to the number of GFP molecules. The absorbance is expressed in dimensionless units, whereas fluorescence intensities have specific relative fluorescence units (RFU). In this section we describe how, by means of the measurement models previously described, we derive promoter activities and protein concentrations from the absorbance and fluorescence data (Berthoumieux et al., 2013b;[START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF].

Reconstruction of promoter activities

The corrected absorbance and fluorescence data are used to compute promoter activities (synthesis rates) and protein concentrations, following the measurement models in [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]. From Eqs. 2.18-2.19 it follows that

f (t) = d dt r(t) + (µ(t) + γ r ) r(t).
(2.25)

The growth rate µ(t) can be estimated from the absorbance, that is,

µ(t) = d dt A(t) 1 A(t) = d ln A(t) dt .
(2.26)

The time-varying GFP concentration in the bacterial population, r(t), can also be estimated from the absorbance and fluorescence, making the usual assumptions that the fluorescence is proportional to the number of GFP molecules and the absorbance proportional to the biomass:

r(t) ∼ I(t) A(t) .
(2.27)

We arbitrarily set the proportionality constant in Eq. 2.16 to 1, thus expressing the reporter protein concentration in units RFU (and the synthesis rate in units RFU min -1 ). Substituting the expressions for r(t) and µ(t) into Eq. 2.25 and after some basic computations (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF] we obtain:

f (t) = d dt I(t) A(t) + γ r I(t) A(t) .
(2.28)

The definition is equivalent to other definitions in the literature [START_REF] Ronen | Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics[END_REF] when µ(t) >> γ r . The expression is evaluated using estimates of A(t), I(t), and dI(t)/dt obtained by means of cubic smoothing splines (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF].

Reconstruction of protein concentrations

In order to reconstruct the concentration of a protein of interest, the same measurement models are used, in particular Eq. 2.17. The term kp g(t) was seen to be proportional to f (t), following Eq. 2.20. We arbitrarily set the proportionality constant in Eq. 2.20

to 1, and we rewrite:

d dt p(t) = f (t) -(µ(t) + γ p ) p(t), p(0) = p 0 , (2.29)
3. Inference of quantitative models of bacterial promoters from timeseries gene expression data This chapter will present the results of the PhD thesis. We will show how the transcriptional response of the genes in the FliA-FlgM module and global regulatory effects have been measured by means of fluorescent reporter genes, in a variety of wild-type and mutant conditions, in different growth media. We will present the mathematical models developed to describe FliA-dependent gene expression. Furthermore, we will illustrate how these data were used to systematically test the information required for the reliable inference of the regulatory interactions and quantitatively predictive models of gene regulation. In a first step, we tested if the use of FliA and FlgM promoter activities, instead of their protein concentrations, allows the expected pattern of regulatory interactions to be inferred, and a quantitative model of the activity of FliA-dependent genes to be identified from the data. In a second step, we introduced global regulatory effects, measured by means of a reporter gene driven by a constitutive promoter.

In a third step we estimated the concentrations of FliA and FlgM from the observed promoter activities and physiologically plausible half-lives of the proteins. The results had been further refined in a fourth step, by taking into account that FliA and FlgM half-lives may vary across conditions, in the range of physiologically valid values.

We also describe in detail the experimental methods used either to produce or validate the biological data on the central module controlling motility in E. coli, along with the experimental conditions, the strains and the inference and modeling frameworks.

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

Results

3.1.1 Monitoring the transcriptional response of the FliA-FlgM module

The more than 60 genes responsible for motility in bacteria are structured in a transcriptional hierarchy of three operon classes, which has been mapped in detail for

Escherichia coli and Salmonella enterica [START_REF] Chevance | Coordinating assembly of a bacterial macromolecular machine[END_REF][START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF][START_REF] Kutsukake | Transcriptional analysis of the flagellar regulon of Salmonella typhimurium[END_REF]Macnab, 1996a). The single class 1 operon flhDC encodes the proteins FlhD and FlhC, that form a heteromultimeric complex activating The transcriptional hierarchy produces a temporally-arranged order of events during the assembly of the flagella and the chemotactic sensing system [START_REF] Chevance | Coordinating assembly of a bacterial macromolecular machine[END_REF][START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF][START_REF] Kutsukake | Transcriptional analysis of the flagellar regulon of Salmonella typhimurium[END_REF]Macnab, 1996a). On the highest level of the hierarchy, the transcription of the flagellar master regulator responds to a variety of intracellular signals [START_REF] Girgis | A comprehensive genetic characterization of bacterial motility[END_REF][START_REF] Pesavento | Inverse regulatory coordination of motility and curli-mediated adhesion in Escherichia coli[END_REF]. For instance, the expression of the flhDC operon is repressed when the bacteria are grown on minimal medium with glucose [START_REF] Adler | The effect of environmental conditions on the motility of escherichia coli[END_REF]. When glucose is depleted from the environment, however, the signalling molecule cyclic AMP (cAMP) accumulates in the cell, which induces flhDC transcription through the intermediary of the cAMP receptor protein Crp [START_REF] Zhao | Adaptation in bacterial flagellar and motility systems: from regulon members to "foraging"-like behavior in E. coli[END_REF]. In the presence of FlhDC, the class 2 operons, and thus the genes encoding the hook basal-body (HBB) structure as well as FliA and FlgM, are actively transcribed. FlgM sequesters FliA and prevents it from transcribing the observation that most FlgM in the cell derives from the FliA-dependent promoter [START_REF] Chevance | Coordinating assembly of a bacterial macromolecular machine[END_REF][START_REF] Gillen | Transcription from two promoters and autoregulation contribute to the control of expression of the salmonella typhimurium flagellar regulatory gene flgm[END_REF]. In the analysis that follows, we therefore neglected flgM transcription from the class 2 promoter.

As illustrated in Figure 3.2, the primary absorbance and fluorescence signals can be transformed into promoter activities using kinetic models of gene expression (Sec-

Results

tion 3.2.2). More precisely, the reporter gene data allow one to deduce protein synthesis rates (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF][START_REF] Ronen | Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics[END_REF]. Under certain conditions, as explained in detail in Chapter 2, the latter are proportional to mRNA concentrations and promoter activities and thus reflect the transcriptional activity of the gene. Following established terminology, we will refer to the measured protein synthesis rates as promoter activities, or more generally, activities of genes. Activity of the tar promoter (•, blue), computed from the primary data as described in Section 3.2.2 and in Chapter 2. The solid line corresponds to the mean of 6 replicate absorbance measurements and the shaded region to the mean of the promoter activities ± twice the standard error of the mean.

In each of the experimental conditions, we have acquired 5 to 8 replicate measurements, which allows for an estimation of the uncertainty in the derived promoter activities. Figure 3.3 shows the results for the five conditions considered here: (i) ∆rpoS strain grown in M9 (∆rpoS -M9), (ii) ∆cpxR strain grown in M9 (∆cpxR-M9), (iii) ∆csgD strain grown in M9 (∆csgD-M9),(iv) ∆csgD strain grown in LB (∆csgD-LB), and (v) wild-type strain grown in LB (WT-LB). As expected (Adler and Templeton,

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

1967), the fluorescence signals in the wild-type strain grown in glucose were mostly not distinguishable from the background fluorescence and therefore this condition was not further considered. In one condition (WT-LB), the activities measured by means of reporter genes were validated using RT-qPCR (Section 3.2.6).

The measured activity profiles in Figure 3.3 show some common features, such as a transient activity peak of the genes during exponential growth, followed by stabilization at a low level after growth arrest. The induction of the individual promoters has a distinct temporal order, corresponding to the level of the promoters in the transcriptional hierarchy [START_REF] Kalir | Using a quantitative blueprint to reprogram the dynamics of the flagella gene network[END_REF]): fliA, flgM, tar. There are also clearly visible differences between the profiles across the conditions though. In M9 medium with glucose the motility genes in the mutant strains are transcribed right from the start, whereas in LB induction occurs only after a number of generations, consistent with previous reports [START_REF] Adler | The effect of environmental conditions on the motility of escherichia coli[END_REF][START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF]. Moreover, the strength of induction and the duration of the activity peak varies from one condition to the other. For instance, the maximal activity of tar varies 10-fold between the WT-LB and ∆csgD-LB conditions.

Identification of gene regulation functions from promoter activities

The circuit in Figure 3.1 has been well-studied over several decades and its regulatory structure is well-known [START_REF] Keseler | EcoCyc: fusing model organism databases with systems biology[END_REF]. This therefore provides an excellent test case for investigating what kinds of information are needed for the reliable inference of regulatory interactions and quantitative regulation functions from gene expression data. In a first step, we tested if we could account for measured time-varying promoter activities while ignoring the distinction between mRNA and protein concentrations as well as the activity of the gene expression machinery and other global physiological effects, as is usually the case.

We expect FliA to be an activator and FlgM an inhibitor of target genes like tar and flgM. In order to check if this regulatory pattern is consistent with the reporter The promoter activities of fliA (green), flgM (red), and tar (blue) measured in all experimental conditions considered in this study: ∆rpoS strain grown in M9 (∆rpoS -M9), ∆cpxR strain grown in M9 (∆cpxR-M9), ∆csgD strain grown in M9 (∆csgD-M9), ∆csgD strain grown in LB (∆csgD-LB), and wild-type strain grown in LB (WT-LB). The promoter activities have been derived from the primary data as illustrated in Figure 3.2.

gene data, we used minimal sign pattern analysis (Porreca et al., 2010a). This approach exploits time-series data to invalidate patterns of regulatory interactions, based on the assumption that the activity of a gene is a monotonic function of its regulators.

The remaining patterns of regulatory interactions are subsumed by so-called minimal patterns. These patterns are minimal in the sense that removing any of the regulators results in an inconsistency with the data, while adding other regulators preserves consistency (see Section 3.2.3 for details on the method).

We applied minimal sign pattern analysis to the reporter gene data in Figure 3.3.

In particular, we tested whether the expected regulatory pattern is conserved when replacing the concentrations of FliA and FlgM by the measured promoter activities. In order to check the robustness of the minimal patterns thus obtained, we verified that no regulatory patterns were dismissed because of a single measurement in the time-

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

series. We found that, both for the tar and the flgM promoter, the expected regulation by FliA and FlgM is not consistent with the data (Figure 3.4). Intuitively, this can be explained by the fact that, over some interval of time in the condition ∆rpoS , a decrease of the promoter activity of fliA and an increase of the promoter activity of flgM coincide with an increase of the activity of the target genes.

FliA FlgM

Tar FliA FlgM FlgM Figure 3.4: Minimal sign patterns for the regulation of tar and flgM when replacing protein concentrations by promoter activities. For both the regulation of tar gene and flgM gene, the expected sign pattern (fliA, flgM, tar ) = (1,-1,0) is found to be inconsistent with the data. The invalidation of the expected sign pattern is due to the fact that, in ∆rpoS, a decrease of the promoter activity of fliA and an increase of the promoter activity of flgM corresponds to an increase of the promoter activity of tar (and flgM ). The minimal sign pattern identified for tar gene is (0,0,1), meaning that tar is necessary for its own regulation. The minimal sign pattern identified for flgM gene is (0,1,0), meaning that flgM is necessary for its own regulation. Black arcs represent the minimal consistent sign patterns. Every consistent pattern can be obtained from one of these minimal sign patterns by turning some gray arcs into black arcs with either a line (inhibition) or an arrow (activation) end.

Despite this structural problem, we also tested to which extent it is possible to quantitatively predict the activities of tar and flgM from the activities of their regulators.

To this end, we developed a mechanistic model of the regulation of these promoters by FliA and FlgM. The model takes into account the titration of FliA by FlgM and the activation of transcription by (free) FliA. We made a quasi-equilibrium assumption for FliA-FlgM association and dissociation, justified by the fast time-scale on which these reactions occur in comparison with transcription and translation processes [START_REF] Buchler | Molecular titration and ultrasensitivity in regulatory networks[END_REF][START_REF] Bundschuh | Fluctuations and slow variables in genetic networks[END_REF]. Moreover, we chose a Hill function to describe promoter activation and included a basal synthesis rate. The resulting model

Results

is:

f (t) = k 0 + k 1 p A,f ree (t) n θ n + p A,f ree (t) n , (3.1) p A,f ree (t) = 1 2 -(K + p M (t) -p A (t)) + (K + p M (t) -p A (t)) 2 + 4 K p A (t) , (3.2)
where f (t) is the time-varying promoter activity, p A,f ree (t) is the concentration of free FliA, θ is a threshold constant for promoter activation, k 0 and k 1 are the basal and maximal synthesis rates, respectively, and n is a Hill constant. The concentration of free FliA is computed from the concentrations p A (t) and p M (t) of total FliA and FlgM, respectively, and the FliA-FlgM dissociation constant K. All variables and parameters are non-negative and n ≥ 1. The concentration variables, as well as θ and K, have the units RFU, while the promoter activity and the rate constants have units RFU min -1 . The derivation of the model is described in detail in Section 3. in the regulation model from the data obtained in all five conditions, using a hybrid genetic algorithm that was shown to give good results for nonlinear models in systems biology [START_REF] Rodriguez-Fernandez | Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems[END_REF]. The algorithm minimizes the mean-square error between the observed promoter activities and the predictions of the model of Eqs. 3.1-3.2, while taking into account differences in absolute promoter activity across conditions as well as the time-varying size of confidence intervals (Section 3.2.2).

The predictions of the identified regulation function for tar are shown in Figure 3.5.

While the fit with the experimental data is quite good for the ∆csgD-LB condition and acceptable for the WT-LB condition, the model is not able to account for the peak in tar activity in the M9 conditions. The model either predicts no peak or a peak occurring more than an hour before it is observed. When analyzing the estimated parameter values, we observe that the cooperativity parameter n equals 1 and that the This means that the the regulation function of the tar promoter is essentially a linear transformation of fliA activity. The parameters were estimated using a multistart global optimization algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 33.6 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) = (7.6, 853, 1, 662, 14615).

The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

We repeated the above analysis for the flgM promoter, setting the parameter that is not promoter-specific, the FliA-FlgM dissociation constant K, to the value estimated from the tar data. Since the fluorescence signal emitted by the strain carrying the pflgM -gfp reporter plasmid is very close to the background levels, and thus unreliable, we eliminate the condition WT-LB. The results are shown in Figure 3.6 and are qualitatively similar to results obtained for tar.

In conclusion, replacing protein concentrations by promoter activities in the FliA-FlgM module is insufficient for obtaining reliable models of the promoter activities, either structurally or quantitatively. The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 23 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) = (9, 582, 1, 221, 7307). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Identification of gene regulation functions from promoter activities including global physiological effects

A possible explanation for the difficulty to identify quantitative regulation functions from information on promoter activities only may be that, in addition to transcription regulators and other specific regulators, the activity of the gene expression machinery also affects gene expression [START_REF] Bremer | Modulation of chemical composition and other parameters of the cell by growth rate[END_REF][START_REF] Klumpp | Growth rate-dependent global effects on gene expression in bacteria[END_REF][START_REF] Maloe | Biological Regulation and Development, chapter Regulation of the protein synthesizing machinery, ribozomes, tRNA, factors and so on[END_REF]. Contrary to FliA and FlgM, which affect specific genes, all motility genes are affected by the activity of the gene expression machinery and other global physiological effects. Figure 3.7 shows the network structure of the FliA-FlgM module when such global physiological effects are taken into account.

The activity of the gene expression machinery includes the abundance and activity of RNA polymerase and ribosome, as well as pools of metabolic precursors, and is therefore difficult to quantify in a direct way. This has motivated the use of the growth rate and the activity of constitutive genes, whose expression is in principle not controlled by any specific regulators, as indirect read-outs of the global physiological state of the cell (Berthoumieux et al., 2013b;[START_REF] Gerosa | Dissecting specific and global transcriptional regulation of bacterial gene expression[END_REF][START_REF] Klumpp | Growth rate-dependent global effects on gene expression in bacteria[END_REF]. In this study, following (Berthoumieux et al., 2013b), we use the activity of the phage λ promoter pRM, which is constitutive in non-infected E. coli cells, as a quantitative measure of the activity of the gene expression machinery, and the global physiological state more generally. In Figure 3.8 the time-varying activity of the constitutively-expressed reporter gene is shown, together with the activity of tar.

Does the inclusion of global physiological effects enable the identification of quantitatively predictive gene regulation functions?

In order to answer this question, we again applied minimal sign pattern analysis to the reporter gene data, this time includ- ing the activity of the constitutive phage promoter as a proxy for the activity of the gene expression machinery. As in the previous section, the FliA and FlgM concentrations were replaced by the activities of their genes. Whereas the expected pattern of regulatory interactions (activation of the promoter by the gene expression machinery and FliA, inhibition by FlgM) was consistent with the data for tar, the analysis again ruled out this pattern for flgM (Figure 3.9). This means that, even when including global physiological effects in the analysis, the regulatory structure cannot generally be recovered.

Ignoring the fact that the correct structure could not be recovered for flgM regulation, we also checked if the proposed extension improves the capability of the regulation function for FliA-controlled promoters to quantitatively account for the time-varying data. To this end, we multiplied Eq. 3.1 with f const (t), the measured activity of a constitutive promoter:

f (t) = f const (t) k 0 + k 1 p A,f ree (t) n θ n + p A,f ree (t) n , (3.3)
The fits shown in Figure 3.10, obtained with the parameter estimation approach The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 30 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) = (0.1, 16, 1.04, 662, 14615). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure .   functions are still unsatisfactory. As explained in Chapter 1, replacing protein concentrations by promoter activities may not be appropriate, due to the fact that the half-lives of proteins are usually much longer than the half-lives of mRNA, causing the temporal decorrelation of protein concentrations and promoter activities. We therefore investigated how information on protein concentrations can be integrated into the inference process.

Identification of gene regulation functions from estimates of protein concentrations

It is straightforward to provide an estimate of the GFP concentration from the fluorescence and absorbance data, as explained in Section 3.2.2. The results are shown in

Time (min)

Figure 3.12: Estimates of FliA concentrations from reporter gene data. Concentrations of FliA (dashed line) computed from the fliA promoter activity (solid line) in all experimental conditions considered in this study. The fliA activities are the same as shown in Figure 3.3. The dark green line represents the concentration of the reporter protein, while the light green line represents the reconstructed concentration for the measured half-live of 30 min. Promoter activity has been normalized with respect to the maximum of the upper limit of its confidence interval in each condition. All protein concentrations have been normalized with respect to the maximum of the upper limit of the confidence interval of the reporter concentration in each condition. The shaded region corresponds to the mean of the promoter activities ± twice the standard error of the mean. Similar estimates of FlgM concentrations can be found in Figure 3.13.

FlgM in Salmonella enterica wild-type strains growing in LB were found to be 30 min and 18 min, respectively [START_REF] Aldridge | The flagellar-specific transcription factor, σ 2 8, is the type iii secretion chaperone for the flagellar-specific anti-σ 2 8 factor flgm[END_REF]. These half-lives are much shorter than those commonly found for proteins in E. coli. This can be explained by the fact that, in addition to being physically degraded, FlgM is secreted from the cell. Moreover, FliA is subject to active degradation by Lon (Figure 3.1).

How can we exploit this information to reconstruct the protein concentration from the promoter activity? As shown in (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF] and Section 3.2.2, if the half-live of the protein of interest is known, then an estimate of its concentration can be reconstructed from the observed promoter activity using a simple kinetic model integrating the effects of protein synthesis and degradation as well as growth dilution of the protein. Figure 3.12 shows the result that is obtained for the FliA concentration, using The plots correspond to the five experimental conditions considered (∆rpoS -M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB, and WT-LB) as well as the pooling of the data sets from all five conditions. The dot in the center of each region in the plots corresponds to a tested combination of half-lives of FliA and FlgM, and thus to specific protein concentration profiles computed from the kinetic model of gene expression (Section 3.2.4). The minimal regulatory patterns were obtained by applying the minimal sign pattern algorithm (Porreca et al., 2010a). The color codes represent the different categories of minimal signal patterns inferred. A region is colored green if the expected regulatory patterns is among the minimal sign patterns returned by the algorithm, and yellow if it is compatible with the returned sign patterns. A region is colored red if none of the returned sign patterns is consistent with the data. Two examples of inconsistent sign patterns are shown. The values of the half-lives are represented in log. Similarly to Figure 3.14, the plots correspond to the four of the experimental conditions considered (∆rpoS -M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB) as well as the pooling of the data sets from all five conditions. The condition WT-LB was not used in the analysis of the regulation of the flgM promoter. The dot in the center of each region in the plots corresponds to a tested combination of half-lives of FliA and FlgM, and thus to specific protein concentration profiles computed from the kinetic model of gene expression (Section 3.2.4). The minimal regulatory patterns were obtained by applying the minimal sign pattern algorithm (Porreca et al., 2010a). The color codes represent the different categories of minimal signal patterns inferred. A region is colored green if the expected regulatory patterns is among the minimal sign patterns returned by the algorithm, and yellow if it is compatible with the returned sign patterns. A region is colored red if none of the returned sign patterns is consistent with the data. Two examples of inconsistent sign patterns are shown. The values of the half-lives are represented in log. 3, 4.6, 2.4, 3030, 223750)) and two other fits for comparable half-lives (dashed lines, .2, 4.7, 2.2, 4535, 222800) and 3, 4.6, 2.4, 2800, 162000) ). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure .   sections, at high concentrations of FliA the tar promoter is maximally expressed, as expected.

Results

FliA FlgM

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

Q = 24.1, (k 0 , k 1 , n, θ, K) = (0
Q = 25, (k 0 , k 1 , n, θ, K) = (0.
The above analysis ignores a particularity of the FliA-FlgM module, namely that the half-lives vary across growth conditions. Generally speaking, in environmental conditions favoring a larger number of flagella, and thus completed HBB structures, the secretion rate of FlgM is higher and therefore the apparent half-life shorter. In mutant strains without HBB structures, and thus no protein secretion, the FlgM half-life is 3 h (Karlinsey et al., 2000a), while in some conditions half-lives up to 7 min were measured 3, 6.1, 2, 3170, 223750)) and two other fits for comparable half-lives (dashed lines, Q = 26.3, (k 0 , k 1 , n, θ, K) = (0.4, 5.9, 2. 3,2358,279500) and 3, 6, 2.1, 2760, 162000) ). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure. [START_REF] Karlinsey | Flk couples flgm translation to flagellar ring assembly in salmonella typhimurium[END_REF], see [START_REF] Aldridge | The flagellar-specific transcription factor, σ 2 8, is the type iii secretion chaperone for the flagellar-specific anti-σ 2 8 factor flgm[END_REF] for intermediate values). The half-life of FliA, the flagellar sigma factor, is also variable. FliA is subject to active degradation by the Lon protease, but stabilized when bound to FlgM (Figure 3.1). This makes its apparent half-life dependent on the concentration of its anti-sigma factor [START_REF] Barembruch | Cellular levels and activity of the flagellar sigma factor flia of escherichia coli are controlled by flgm-modulated proteolysis[END_REF]. The measured half-life of FliA in mutant strains without HBB structures, and thus with maximal protection of by FlgM, is 2 h [START_REF] Barembruch | Cellular levels and activity of the flagellar sigma factor flia of escherichia coli are controlled by flgm-modulated proteolysis[END_REF]. In wild-type strains exponentially growing in LB medium, this value may decrease down to 30 min [START_REF] Aldridge | The flagellar-specific transcription factor, σ 2 8, is the type iii secretion chaperone for the flagellar-specific anti-σ 2 8 factor flgm[END_REF]. In summary, the half-lives of

Q = 27.1, (k 0 , k 1 , n, θ, K) = (0.

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

both FliA and FlgM are not identical across all growth conditions considered. While we can give upper and lower bounds on the half-lives, we do not know the exact value in most conditions.

This specific property of the FliA-FlgM module suggests a final extension of the analysis to improve the inference results. We allowed the FliA and FlgM half-lives to vary between physiologically possible bounds in each of the conditions and estimated not only the parameters of the regulation functions, but also the half-lives. In order to reduce the computational complexity of this procedure, we discretized the space of possible half-lives, selecting 27 values each for FliA and FlgM, and we precomputed the protein concentration profiles for each half-life in each of the experimental conditions.

The resulting time-course patterns were used for the same analysis as above. This is also intuitively expected, as FlgM is actively exported in the exponential growth phase considered. The best fit finds a cooperativity parameter equal to 1.8 (or 1.09 in the case of the second fit considered). A priori, positive cooperativity is not expected to occur, since FliA does not form a dimer and has only a single binding site in the promoter region. [START_REF] Buchler | Molecular titration and ultrasensitivity in regulatory networks[END_REF] have shown, however, that the titration of a transcription factor may indeed lead to positive cooperativity [START_REF] Buchler | Molecular titration and ultrasensitivity in regulatory networks[END_REF]. The conditions they indicate in their analysis (p M > K and p M ≈ p A ) are satisfied here. Like for the fit with the measured half-lives in Figure 3.16, the activity of tar varies from 0 to its maximal possible value (

Time (min)

Figure 3.18: Fits of regulation function of tar to reporter gene data when reconstructing protein concentrations from the reporter gene data for physiologically realistic half-lives and including global physiological effects. As in Figure 3.16, but the half-lives have now also been estimated from the data, within a physiologically plausible range. Two example fits are shown, namely the best fit for estimated halflives of FliA and FlgM (solid line, Q = 20.9, (k 0 , k 1 , n, θ, K) = (0.2, 5.1, 1.8, 3145, 17204)) and another example of a high-ranking fit (dashed line, Q = 21.09, (k 0 , k 1 , n, θ, K) = (0.12, 8.5, 1.09, 24566, 88350)). In the case of the best fit, the half-lives of FliA are equal to (50,24,24,35,45) min in the (∆rpoS, ∆cpxR, ∆csgD-M9, ∆csgD-LB, WT-LB) conditions, respectively, while the half-lives of FlgM are equal to (27,18,24,18,18) min. In the case of the second fit, the half-lives of FliA are equal to (60, 30, 24, 60, 30) min and the half-lives of FlgM are equal to (9, 11, 24, 45, 7) min in the above experimental conditions, respectively. The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

The improvement in the fit is less evident in the case of flgM promoter (results are reported in Figure 3.19). One possible explanation is that the fit (in Figure 3.11) obtained using the promoter activities of fliA and flgM was already quite good, due

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

to the fact that the activity of flgM is among the regulators. However, the parameter values obtained when using protein concentrations of FliA and FlgM are similar to those obtained for tar promoter analysis (Figure D.2-F in Appendix D). However, the K parameter value is found to be approximately 10 times the maximum of FlgM concentration. This may be a consequence of the fact that the best fit is achieved for combination of half-lives and K values identified for tar regulation. Similarly to Figure 3.18 the half-lives have been estimated from the data, within a physiologically plausible range. Two examples of fits are shown, namely the best fit for estimated half-lives of FliA and FlgM (solid line, Q = 25, (k 0 , k 1 , n, θ, K) = (0.45, 5.9, 2.4, 2930, 222850)) and another example of a high-ranking fit (dashed line, .4, 6.1, 2.3, 3030, 279550)). In the case of the best fit obtained, the half-lives of FliA are equal to (24 min, 35 min, 27 min, 50 min) in the (∆rpoS, ∆cpxR, ∆csgD-M9, ∆csgD-LB) conditions, respectively. The half-lives of FlgM are equal to (27 min, 11 min, 11 min, 9 min) in the (∆rpoS, ∆cpxR, ∆csgD-M9, ∆csgD-LB) conditions, respectively. In the case of the second fit, the half-lives of FliA are equal to (24 min, 35 min, 27 min, 27 h) and the half-lives of FlgM are equal to (18 min, 13 min, 20 min, 27 min) in the above precised experimental conditions, respectively. The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Q = 25.4, (k 0 , k 1 , n, θ, K) = (0
The reconstruction of protein concentrations from transcription data results in much

Results

better inference results for the FliA-FlgM module. The computation of the protein concentrations requires a simple kinetic model, accounting for protein synthesis and degradation, as well as estimates of the protein half-lives. While this increases the complexity of the data analysis procedures, it reflects the actual dynamics of gene expression and is thus critical for exploiting time-series measurements. Moreover, the availability of information on protein half-lives may not be constraining in practice, since even rough half-live estimates from the literature were seen to preserve the expected interaction pattern and provide a significant improvement of the ability of the models to quantitatively describe the time-varying promoter activity.

Determination of conditions in which protein half-lives and global physiological effects are important

The importance of accounting for global physiological effects and protein half-lives was demonstrated above for the regulation of the expression of tar. The same analysis was repeated for the regulation of the flgM promoter. We found that, for this promoter, the improvement in the fit to the experimental data obtained by including global physiological effects and protein kinetics is much less pronounced than for tar. One possible explanation is that the flgM activity profile happens to be already well explained using the promoter activities of fliA and flgM as proxies for the corresponding protein concentrations (Figure 3.6), thus leaving little space for improvement. In addition, from a mathematical viewpoint, we notice that using the promoter activity of flgM for the fitting of the same quantity may render the regression problem degenerate. Still, these results raise a more general question: When is it important to take into account protein half-lives and global physiological effects?

To answer this question we performed an in-silico analysis where the regulation model of Eqs. 3.2-3.3 is simulated for different protein half-lives and varying strength of the global physiological contribution, using the pfliA, pflgM, and pRM activity profiles reported in Figures 3.3 and 3.8. Identification is then attempted from the simulated data with models ignoring protein half-lives and global physiology. This enables us to quantify the relevance of the analysis in the previous sections for a variety of realistic scenarios, starting from experimentally measured activities of bacterial promoters.

To evaluate the importance of protein half-lives, we simulated FliA and FlgM concentration profiles for half-lives ranging between 7 minutes and 16 hours. The other

f (t) = α • f const (t) -fconst + fconst • k 0 + k 1 p A,f ree (t) n θ n + p A,f ree (t) n , (3.4)
with p A,f ree (t) computed from the FliA and FlgM concentration profiles according to Eq. 3.2. The upper bound of 1.25 for α was chosen so as to avoid negative values of the promoter activity f (t).

Identification results using FliA and FlgM concentrations computed for the reference half-lives of 30 min and 18 min, respectively, but ignoring global physiological effects

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

are reported in Figure 3.20B. It is clear that the misfit of the tar promoter activity data increases with the strength of the ignored physiological effects. In particular, with the experimentally observed pRM activity (α = 1), the discrepancy between the data and the best model fit is quite significant. This is in agreement with the results of previous sections, where accounting for global physiological effects turned out to be an important step toward improving the quality of the inference results. Neglecting small variations of global physiological state (α 1) is safer, but ignoring highly varying global physiological effects (α > 1) may have even more severe repercussions on the inference results than those observed in the previous section.

In summary, the simulation study shows that, as expected, the importance of accounting for protein kinetics and global physiological effects depends on the strength of these effects, although the structure of the system itself may also play a role, as illustrated by the different dependency of the quality of the fit on FliA and FlgM concentrations (Figure 3.20A). As a general rule, ignoring significant fluctuations of the global physiology or large differences between mRNA and protein half-lives is very likely to result in modelling bias and hence poor inference results. Interestingly, in the previous sections a substantial improvement of the fit of a quantitative regulation function to tar activity was already obtained when taking into account concentrations of short-lived proteins and moderately-variable global physiological effects. In the light of the analysis of this section, the contribution of our approach becomes even more fundamental in other systems, bearing in mind that the vast majority of bacterial proteins are much more stable than FliA and FlgM, which are actively degraded and exported from the cell (Figure 3.1).

In the next chapter, we propose further guidelines for experimental design to facilitate the implementation of the approach developed here.

Methods and materials

3.2 Methods and materials

Strains and growth conditions

The E. coli strains we used in this study are all derived from the wild-type strain BW25113. In particular, we used the ∆rpoS , ∆cpxR and ∆csgD deletion strains of BW25113 taken from the Keio collection [START_REF] Baba | Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection[END_REF]. The mutants were reconstructed in our laboratory [START_REF] Dudin | Expression dynamics of rpos/crl-dependent genes in Escherichia coli[END_REF] in order to eliminate the kanamycin resistance gene present in the original deletion strains (Table 3 The wild-type and mutant strains were transformed with low-copy plasmids bearing a gfp reporter gene (Table 3.2). The reporter plasmids for the genes tar, fliA, flgM, and flgA were selected from the plasmid library developed at the Weizmann Institute [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF]. These low-copy pUA66gfp plasmids carry the kanamycin resistance gene and have the origin of replication of the pSC101 plasmid. The promoter regions of the genes of interest control the transcription of the gene encoding the stable GFP-mut2 reporter. The same vector was used to construct a reporter for the constitutive promoter pRM of the phage lambda, by cloning the pRM promoter region contained on the pZE1RMgfp plasmid used in Berthoumieux et al. (2013b) into the pUA66gfp plasmid backbone. Table E.1 in Appendix E lists the primer sequences used for the construction of this pUA66pRM-gfp plasmid using the Gibson Assembly method (Gibson, 2011). The above-mentioned plasmids were transformed into the wild-type and deletion strains of Table 3.1. We verified that the plasmids do not modify the growth of the transformed strains. All strains and plasmids were verified by PCR.

The pRM promoter fused with the gfp reporter gene was also inserted into the chromosome of the BW25113 wild-type strain as reference for the qRT-PCR assays. The WTpRM strain was constructed by using a linear DNA recombination protocole of [START_REF] Sharan | Recombineering: a homologous recombination-based method of genetic engineering[END_REF]. The pRM promoter region along with the gene encoding the GF-Pmut3 reporter were introduced into the intS loci on the chromosome of the BW25113

WT strain, by means of the λ Red system. pRM-gfpmut3 was recovered from the pZE1RMgfp plasmid used in Berthoumieux et al. (2013b).The recombineering protocols use the bacteriophage λ Red system that includes the phage recombination genes gam, bet and exo. The protein coded by gam, Gam, prevents E. coli nuclease from degrading linear DNA fragments [START_REF] Karu | The gamma protein specified by bacteriophage lambda. Structure and inhibitory activity for the RecBC enzyme of Escherichia coli[END_REF][START_REF] Murphy | Lambda gam protein inhibits the helicase and chi-stimulated recombination activities of Escherichia coli recbcd enzyme[END_REF] thus allowing preservation of transformed linear DNA in vivo. The bet gene product, Beta, is a ssDNA binding protein that promotes annealing of two complementary DNA molecules [START_REF] Karakousis | The beta protein of phage lambda binds preferentially to an intermediate in dna renaturation[END_REF], and the exo gene product, Exo, has a 5 to 3 dsDNA exonuclease activity [START_REF] Cassuto | Role of exonuclease and beta protein of phage lambda in genetic recombination. v. recombination of lambda dna in vitro[END_REF]. Working together these latter two proteins insert linear DNA at the desired target, creating genetic recombinants.

For all experiments, the strains were recovered from glycerol stock and grown overnight (16 h) at 37 • C in LB rich medium or M9 minimal medium [START_REF] Miller | Experiments in Molecular Genetics[END_REF] supplemented with 0.3% glucose and mineral trace elements. For the preculture of strains containing plasmids, kanamycin (50 µg/ml) was added. The overnight cultures were diluted (10-to 100-fold) into a 96-well microplate, so as to obtain an adjusted initial OD 600 of 0.2. The wells of the microplate contain 150 µl of the above medium, to which was added 1.2% of the buffering agent HEPES (4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid) for maintaining a constant external pH. The wells were covered with 60 µl of mineral oil to prevent evaporation. The microplate cultures were then grown for about 16 h at 37 • C, with agitation at regular intervals, in a microplate reader (Fusion Alpha, Perkin-Elmer).

Methods and materials

Experimental monitoring of gene expression in real time and data analysis

The expression of the fluorescent reporter genes in different genetic backgrounds and different growth media was monitored in vivo and in real time. About 150 readings each of absorbance and fluorescence were obtained during a typical experiment using the Perkin-Elmer microplate reader. The absorbance measured at 600 nm quantifies the biomass, while the fluorescence signal emitted at 520 nm, when excited at 485 nm, is proportional to the number of GFP molecules. In order to compute promoter activities and protein concentrations from these data, data analysis procedures were designed and implemented in MATLAB, completing earlier work (Berthoumieux et al., 2013b;[START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF][START_REF] Ronen | Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics[END_REF]. These analysis procedures take into account for the specific half-life of the fluorescent reporter protein and implement procedures for subtracting the autofluorescence background.

Background subtraction

We first corrected the absorbance for the background absorbance of the growth medium.

The corrected absorbance signal A(t) is computed as

A(t) = A u (t) -A b (t), (3.5)
where A u (t) is the primary absorbance signal and A b (t) the absorbance of the growth medium (M9 or LB, depending on the experiment).

The fluorescence signal was corrected for autofluorescence generated by wild-type bacteria carrying the pUA66gfp plasmid without any promoter driving the expression of gfp or no plasmid at all (in practice these two measures of the autofluorescence gave the same result). The autofluorescence depends on the (time-varying) population size.

Since the culture generating the fluorescence signal of interest and the culture generating the autofluorescence signal may not be exactly synchronized, direct subtraction of the autofluorescence background is not always possible. We used a calibration procedure, such that the corrected signal I(t) is defined by

I(t) = I u (t) -s(A(t)), (3.6)
where I u (t) is the primary fluorescence level and s a calibration function, mapping absorbance levels to autofluorescence levels. The calibration function is obtained by

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

fitting a cubic smoothing spline to the autofluorescence generated by bacteria carrying the promoterless pUA66gfp plasmid or no plasmid at all as a function of the absorbance.

Splines have the advantage that they can be evaluated for any absorbance within the observed range and easily extrapolated beyond this range.

Computation of promoter activity and protein concentrations

Following the measurement model in (de [START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF], we describe the expression of the gene of interest and of the reporter protein as follows (Chapter 2):

d dt p(t) = f (t) -(µ(t) + γ p ) p(t), p(0) = p 0 , (3.7)
d dt r(t) = f (t) -(µ(t) + γ r ) r(t), r(0) = r 0 , (3.8)
where p(t) and r(t) are the concentrations of the protein of interest and of the reporter protein, respectively, µ(t) is the time-varying growth rate, and γ p , γ r [min -1 ] are the degradation constants of the protein of interest and the reporter protein, respectively.

Notice that in the case of FlgM, protein degradation includes both physical degradation of the protein and secretion through the cell membrane. The reporter protein concentration r(t) and the promoter activity f (t) can be computed by means of the formulas:

r(t) = I(t) A(t) , (3.9)
f (t) = d dt r(t) + (γ r + µ(t)) r(t) = d dt I(t) A(t) + γ r I(t) A(t) , (3.10)
The reporter concentration is expressed in units RFU and the promoter activity in units RFU min -1 , as is usual for this kind of measurements (see (Berthoumieux et al., 2013b) and Chapter 2). The growth rate is easily estimated from the time-varying absorbance, using the standard relation µ(t) = d ln A(t)/dt.

We used cubic smoothing splines (csaps function in MATLAB) to fit the fluorescence and absorbance data and obtain estimates of A(t), I(t), dA(t)/dt, and dI(t)/dt. The half-life of the GFPmut2 reporter used in this study is 18 h (γ r = 0.0006 ± 0.0001)).

The maturation time of GFPmut2 is short enough (4 min, [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF]) to be computed for a half-life of 2 h from the corrected absorbance (-, grey) and corrected fluorescence measurements using Eq. 3.7. F: Protein concentration of tar (•, blue) computed for a half-life of 18 h from the corrected absorbance (-, grey) and corrected fluorescence measurements using Eq. 3.7.

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

ignored.

In order to reconstruct the concentration of a protein of interest, we again use Eq. 3.7.

The promoter activity, f (t), is proportional to the synthesis rate, as explained in detail in Chapter 2. When the degradation constant is known, we can compute the protein concentration by numerical integration, starting from the initial concentration p 0 .

This initial concentration is obtained from the reporter gene data, by realizing that the bacterial cells at the beginning of the experiment are rediluted cells from a preculture grown in the same medium. In particular, assuming that gene expression in the preculture is at steady-state, it follows from Eq. 3.7 and Eq. 3.8 that p(0) = p(T ) = µ(T ) + γ r µ(T ) + γ p r(T ), (3.11) where µ(T ) is the growth rate of the preculture at the time of redilution (at the time T ), p(T ) and r(T ) are the corresponding concentrations of the protein of interest and reporter protein, respectively. Usually, the bacteria in the preculture are in stationary phase, so µ(T ) = 0. Eq. 3.7 was solved by numerical integration using the quad function in Matlab.

In the case of the motility network there are two complications that slightly modify this general scheme. First, the half-lives of FliA and FlgM are variable over the time-course of the experiment. During exponential growth, when the motility genes are expressed, FliA and FlgM have short half-lives, due to proteolysis and secretion, respectively.

During stationary phase, at the end of the preculture, this is no longer the case and FliA and FlgM have larger half-lives. As a consequence, when computing the initial protein concentrations from the reporter concentrations at time T , we need to take protein degradation constants γ p ' corresponding to these larger half-lives. Second, in some experimental conditions, notably in rich medium like LB, the activity of the fliA, flgM, and tar promoter is negligible in the first few hours of the experiment [START_REF] Kalir | Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria[END_REF]. As a consequence, the fluorescence intensity in the corresponding reporter strains is indistinguishable from the background fluorescence. We assume the promoter activity of the genes to be 0 in this case and back-extrapolate the observed promoter activities at earlier times towards 0. For each of the derived quantities r(t), f (t), and p(t), confidence intervals (defined as ±2 standard errors of the mean) were computed from 6-7 experimental replicates.

Computation of minimal consistent sign patterns of regulatory interactions

In this section we adopt a notation, where vector x = (x 1 , . . . , x n ) indicates concentration of regulators (e.g., FliA and FlgM), but may comprise in addition regulatory effects such a global cell physiology, depending on the context. What follows applies identically to all target genes of interest, i.e. tar and flgM. Let f (x) be the promoter activity of one gene of interest.

We use the approach introduced in (( Porreca et al., 2010a)), which exploits time-series data of protein concentrations and promoter activities (protein synthesis rates) to infer patterns of regulatory interactions. The method relies on two assumptions:

1. f (x) is monotonic in every x j , with j = 1, . . . , n; Assumption 1 of the method reflects the hypothesis that a regulator (e.g., a transcription factor, but also the gene expression machinery) cannot operate both as a repressor and as an activator of a specific target gene (see Porreca et al. (2010a) and references therein), while it is allowed to operate as a repressor for one gene and as an activator for another gene. This corresponds to assuming that the activity of a gene is a monotone nondecreasing function of activators and a monotone nonincreasing function of repressors. Any such regulatory pattern can be encoded in terms of a sign pattern, The observed promoter activtity of flgM (dashed line, red) and its etrapolation (*, red). B : The observed promoter activtity of fliA (dashed line, green) and its etrapolation (*, green). C : The effect on the computation of the protein concentration of FlgM when taking into account the extrapolation of its promoter activity in A and the initial condition. Protein concentration was computed for an initial half-life of 3 h and a short half life of 18 min. D: The effect on the computation of the protein concentration of FliA when taking into account the extrapolation of its promoter activity in B and the initial condition. Protein concentration was computed for an initial half-life of 2 h and a short half life of 30 min. The hatched regions correspond to the regions where the activity of promoters was extrapolated, i.e. [0, 410 min] for flgM and [0, 340 min] for fliA. The promoter activities have been derived from the primary data as illustrated in Figure 3.2. The shaded regions correspond to the mean of the promoter activities and protein concentrations, respectively, ± twice the standard error of the mean. The absorbance is drawn in solid, grey lines.

i.e., a vector containing one entry per regulator, taking value +1 for activators, -1

for repressors, and 0 for factors that do not affect the expression of the gene under consideration. We may thus define the sign pattern π = (π 1 , . . . π n ) of f by posing The rationale of the procedure for eliminating hypotheses from the set of all candidate sign patterns is the following (Porreca et al., 2010a). Given any two concentration vectors x and x , the implication

π j = 1 if f is increasing in x j , π j = -1 if f is decreasing in x j ,
π j (x j -x j ) ≥ 0, j = 1, . . . , n ⇒ f (x ) ≥ f (x )
is satisfied by the very definition of the sign pattern π of f . Therefore, for a hypothetical sign pattern π and perfect measurements ( k = e k = 0 for all k), any two data points (x , f ) and (x , f ) that falsify the implication allow one to conclude that π is not the sign pattern of f . In particular, if f < f , the sign pattern π defined by πj = 1 if

x j > x j , πj = -1 if x j < x j , and πj = 0 otherwise, is inconsistent with the data. In addition, any subpattern of π, i.e. a pattern π whose nonzero entries are equal to the corresponding entries of π (denoted with π π), is also inconsistent with the data, since the implication above is still violated.

For instance, in the network module considered in this paper, the assumption that both FlgM and FliA activate tar can be rejected if two measurement times are found such that, for increasing concentrations of FlgM and FliA, the promoter activity of tar is decreased. The algorithm makes the above verifications in a computationally efficient way and returns, for every target gene, a set of minimal sign patterns. The minimal sign patterns are regulatory patterns consistent with the data, having the properties that removal of any interaction results in an inconsistent pattern, whereas addition of a regulator (activator or repressor) preserves the consistency. This test is easily robustified to account for measurement uncertainties, see Figure 3.23 for a graphical example on a network resembling tar regulation. • For all pairwise different data points (x , f ) and (x , f ) in D:

If f < f -Define π = (π 1 , . . . , πn ) by πj =      1, x > x -1, x < x 0, otherwise , with j = 1, . . . , n -Include π in Π • Return Π
At this stage, a generic pattern π is inconsistent if and only if π π for some π ∈ Π (Porreca et al., 2010a).

Computation of minimal consistent patterns Π * from Π

• Set Π * = ∅

• For c = 0, 1, . . . , n:

-Enumerate all possible patterns π of complexity c

-For every such π:

If π ∈ Π such that π π (π is consistent), and

If π * ∈ Π * such that π * π (π is minimal consistent), then
Include π in Π *

• Return P *

Methods and materials

x 1

x 0 1 ' e 0 1

x 00 1 ' e 00 1

x 2

x 0 2 ' e 0 2

x 00 2 ' e 00 2 x 3

x 0 3 ' e 0 3

x 00 3 ' e 00

1 2 3 3 A B C D E 1 1   0 1   0 1   0 1   0 2   0 2   0 2   1 2    1 3    0 3   0 3   1 3   Figure 3
.23: Computation of inconsistent and minimal consistent sign patterns from data. An example of the method of Section 3.2.3 is shown for the regulation of gene 3 in a hypothetical network with genes 1,2,3. A: From top to bottom, example time profiles and corresponding confidence intervals (thin black lines) for the concentrations of the proteins enconded by genes 1,2,3 and the synthesis rate of gene 3. Only the two data points (x , f 3 ) and (x , f 3 ) are considered in this example. Non-overlapping confidence intervals of f 3 and f 3 (reported next to each other for ease of comparison by the orange and blue shaded regions) imply f 3 < f 3 . Similarly, non-overlapping confidence intervals for x 1 , x 1 and for x 2 , x 2 imply x 1 < x 1 (π 1 = -1) and x 2 > x 2 (π 2 = 1), respectively, while confidence intervals for x 3 and x 3 overlap (π 3 = 0). Whence, π = (-1, 1, 0). If this was the sign pattern of f 3 , then f 3 (x) should increase for x 1 decreasing to x 1 and x 2 increasing to x 2 (x 3 is irrelevant in the hypothesis π3 = 0), whereas the observation says that f 3 < f 3 . Pattern π = (-1, 1, 0) is thus inconsistent with the data. B -E: Regulation patterns for gene 3, corresponding to the consistent sign patterns of f 3 deduced from the inconsistent patterns Π = {π} obtained in A. Circles represent genes; directed arcs represent regulation of the target gene 3 by regulator j, with j = 1, 2, 3; arrow ends represent activation (π j = 1), line ends represent inhibition (π j = -1). Black arrows represent the minimal consistent sign patterns (1, 0, 0) in B, (0, -1, 0) in C, (0, 0, 1) in D, and (0, 0, -1) in E. Every consistent pattern is obtained from one of the cases B -E by turning the corresponding π j = 0 into either π j = -1 or π j = 1.

At this stage, a generic pattern π is consistent if and only if π * π for some π * ∈ Π * (Porreca et al., 2010a).

Methods and materials

with n ≥ 1. Multiplying by maximal synthesis rate k 1 and adding basal (unregulated) synthesis rate k 0 leads to the model we will be using to describe regulation of fliAdependent genes

f (t) = k 0 + k 1 p n A,f ree p n A,f ree + θ n , (3.13)
Note that, in accordance with the expected regulatory pattern, the function

k 0 + k 1 p n A,f ree (p A ,p M ) p n
A,f ree (p A ,p M )+θ n is increasing in p A and decreasing in p M . To verify this, it suffices to show that derivatives with respect to p A and p M are nonnegative and nonpositive, respectively. Ignoring k 0 and k 1 without loss of generality, the derivative of p A,f ree with respect to p A can be written as

1 2 + 1 2 • -(p M -p A + K) + 2K (p M -p A + K) 2 + 4Kp A . This expression is obviously positive if p M -p A + K ≤ 0. If instead p M -p A + K > 0,
note that the expression is still positive if the square of the fraction,

(p M -p A + K) -2K 2 (p M -p A + K) 2 + 4Kp A ,
is smaller than 1. But this is apparent since, under p M -p A + K > 0, the numerator is no bigger than (p M -p A + K) 2 , whereas the denominator is no smaller than the same quantity. Similarly, the derivative of p A,f ree with respect to p M can be written as

- 1 2 + 1 2 • p M -p A + K (p M -p A + K) 2 + 4Kp A . This expression is obviously negative if (p M -p A +K) ≤ 0. If instead (p M -p A +K) > 0,
note that the square root is no smaller than (p M -p A +K), hence the rightmost fraction is no bigger than 1, i.e. the overall expression is again negative.

The additional regulatory effect of the global physiological state of the cell is quantified via further multiplication by a function f const (t):

f (t) = f const (t) k 0 + k 1 p n A,f ree p n A,f ree + θ n , (3.14)
Monotonicity with respect to p A and p M remains unchanged. In addition, f is increasing in f const . In all cases, the model depends on the (nonnegative) parameters k 0 , k 1 , n, θ, K.

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

Parameter estimation

The promoter activity models we have considered in the Section 3.1 have the form

f (t) = f x(t), c
, where c is a vector of unknown parameters and x is a vector of regressors. The specific form of f x(t), c is given for tar in Eq. 3.13 and 3. Given measurements xs (t), f s (t) of x(t), f (t) (averages of 6-7 experimental replicates) at times t ∈ T s along with confidence intervals ( f s (t) ± s (t)) (computed from the same experimental replicates with s equal to twice the standard error of the mean f s ), we estimate c by solving the optimization problem

min c∈C Q(c), Q(c) = s∈ § t∈ s 1 2 s (t) f s (t) -f x(t), c 2 .
The solution is found in MATLAB using the numerical global search function gs with standard settings (interior-point method fmincon for local minimizations). For tar activity, the search is initialized at the values (k 0 , k 1 , n, θ, K) = ( k0 , k1 , n, θ, K) defined as k0 = min{ f s (t) : t ∈ T s , s ∈ S}; k1 = max{ f s (t) -k0 : t ∈ T s , s ∈ S}; n = 1; K = x1 , where the double bar stands for mean over t ∈ T s and s ∈ S; and, in view of Eq.3.12,

θ = 1 2 -( K + x2 -x1 ) + ( K + x2 -x1 ) 2 + 4 K x1 .
The parameter search space C is given by the constraints

k 0 ≥ 0, k 1 ≥ 0, n ∈ [1, 3], θ ∈ [0, 10 × pA,free ], K ∈ [0, 10 × K max ], where K max = max{x s 2 (t) : t ∈ T s , s ∈ S}.
For the estimation of the regulation function of flgM, the condition WT-LB is not available and hence excluded from the computation of Q(c). Moreover, K is fixed for biological consistency, in this case, to the value inferred from the fitting of tar promoter activity.

Validation of reporter gene data using quantitative RT-PCR

We verified the reporter gene measurements by means of qRT-PCR in the WT-LB condition, following a previously described protocol [START_REF] Lee | Absolute and relative QPCR quantification of plasmid copy number in Escherichia coli[END_REF].

Methods and materials

According to Eq. 2.13 in Chapter 2, the ratio of the promoter activities f 1 , f 2 of two genes is proportional to the ratio of the mRNA concentrations m 1 , m 2 , that is, .15) Measuring gene expression by qRT-PCR allows the relative abundances of the mRNA of a target gene to be quantified with respect to the mRNA of a reference gene ( Van-Guilder et al., 2008). This provides a direct way to verify if the relative promoter activities measured with reporter genes are confirmed by another, independent experimental method. We compared the promoter activity of tar, as an example of a motility gene, with the activity of the constitutive pRM promoter. The validation of the ratio f tar /f pRM was carried out by means of the WTpRM strain, a modified BW25113 strain carrying a natural copy of tar and a transcriptional fusion of the pRM promoter with a gfp gene inserted into the intS loci on the chromosome (Section 3.2.1). Quantitative RT-PCR was used to quantify the relative abundances of tar and gfp mRNA, using a standard qPCR protocol [START_REF] Lee | Absolute and relative QPCR quantification of plasmid copy number in Escherichia coli[END_REF].

f 1 (t) f 2 (t) = k p,1 k p,2 m 1 (t) m 2 (t) . ( 3
We took 5 µL samples at 11 time-points from cultures of the WTpRM strain, growing in a microplate under the conditions described in Section 3.2.1. Total mRNA was protected using the RNAprotect Bacteria Reagent kit (Quiagen) and then extracted using the RNeasy mini kit (Quiagen) according to the protocols of the manufacturer.

The RNA samples were then treated using the turbo DNAse (Ambion) to avoid DNA [START_REF] Lee | Absolute and relative QPCR quantification of plasmid copy number in Escherichia coli[END_REF].

The results were analyzed by means of a standard model for computing m 1 (t)/m 2 (t) at the sample time-points t with respect to m 1 (t 0 )/m 2 (t 0 ), the same quantity at a ref-

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

erence time-point t 0 [START_REF] Pfaffl | A new mathematical model for relative quantification in real-time rt-pcr[END_REF]:

q(t) = m 1 (t)/m 1 (t 0 ) m 2 (t)/m 2 (t 0 ) = E ∆C gf p T gf p E ∆C tar T tar , ( 3
(t) = C gf p T (t)-C gf p T (t 0 ), ∆C tar T (t) = C tar T (t)-C tar T (t 0 ).
As our reference time-point, we chose a measurement during exponential growth on glucose. As a consequence, the changes in mRNA abundance are relative to the mRNA abundance in exponential phase. The efficiencies were measured to be 109% for gfp (E gf p = 2.09) and 105% for tar (E tar = 2.05).

From Eq. 3.15 it follows that

q(t) = f 1 (t)/f 1 (t 0 ) f 2 (t)/f 2 (t 0 ) . (3.17)
The right-hand side of this equation can be computed from the measured promoter activities, as explained in Chapter 2. Figure 3.24 compares the value of q(t) measured by means of reporter genes and qRT-PCR. There is a good qualitative and quantitative correspondence between the two independent methods (qRT-PCR and gene reporter genes) for measuring gene expression.

Conclusion 4.1 Summary of results

In this thesis, we develop methods for inferring genetic regulatory interactions and quantitative gene regulatory functions from gene expression profiles. The validation of the approach proposed in this work is carried out on a well-studied but complex gene regulatory network responsible for the motility of the bacterium E. coli . This chapter summarizes the main contributions of the thesis and puts our results into perspective relative to other methods for network inference.

While many problems persist in current methods for inferring gene regulatory networks from expression data, we believe that the solution should not only be sought in technical improvements of the algorithms themselves, but may come from a better understanding of the precise information on gene expression provided by the experimental data and their integration into appropriate modeling formalisms. The relation between the primary data and physiological quantities like the cellular concentrations of mRNA and protein is usually indirect and obscured by simplifications and assumptions that do not generalize beyond the specific situations for which they were designed. The main regulators of gene expression are proteins. Even though the concentrations of mRNA and proteins are weakly correlated at steady state, this is generally not the case when considering time-varying, dynamical expression data.

Reporter gene systems can yield gene expression profiles in complex in vivo experiments. Such dynamical data are well described by ordinary differential equations. We have therefore adapted an ODE modeling framework to convert (fluorescent) reporter gene data into biologically relevant quantities, such as promoter activities and protein concentrations. This approach is more appropriate for the inference of gene networks in the sense that it explicits the relation between experimental data and physiological quantities by means of mathematical (measurement) models of gene expression. Although an important source of the incomplete and sometimes spurious findings of classic inference algorithms, not many modeling frameworks for inference of gene regulatory networks focus on the distinction between these quantities.

CONCLUSION

Computing protein concentrations by means of the measurement models presented in this thesis supposes that we know the approximate values of the protein half-lives.

With the exception of yeast, genome-wide studies of the stability of individual proteins in microorganisms are rare. It should be noted though that most proteins in E. coli are stable, with half-lives >10 h, so that the decay of protein concentrations is dominated by growth dilution, that is, µ γ p . In other words, in many situations, in order to obtain a reasonable estimate of the effective protein half-life, it will be sufficient to perform the experiments in growth media supporting bacterial growth rates that results in doubling times well below 10 h.

Moreover, many methods for network inference rest on the common and tacit belief that the regulation of gene expression in bacteria is controlled uniquely by transcription factors and other specific regulators. In fact, the complex regulatory activity of the To provide an integrated and straightforward inference approach, our method combines the above described experimental and computational models with a structural and parametric identification algorithm. The algorithm aims at inferring the structure of gene networks from time-series data by exploring the monotonicity properties of the network (or the model) and recovering only structures in good agreement with the data.

For these structures, parameter estimation is performed to find the best prediction of both the (qualitative) model structure and the (quantitative) parameters. Generally, large search spaces reduce the performance, and may even compromise results of inference algorithms, the selection of data-consistent gene network models allows us to focus on precisely analyzing a small set of candidate models. Ultimately, this leads to results that are interpretable and relevant to the initial biological question.

The practical validation of our approach rests on the study of a gold standard biological system, the FliA-FlgM (motility) module in E. coli. Although atypical, this module possesses rich dynamics. The short half-lives of the FliA and FlgM proteins are time-varying and depend on the experimental conditions, inducing the time course of flagella synthesis. We investigated the capacity of our approach to infer from reporter gene data both the regulatory structure and the quantitative regulation function of two uniquely FliA-dependent motility genes (flgM and tar ). When integrating information on the activity of the gene expression machinery and reconstructing protein concentrations from promoter activities, both the structure and the dynamics of the regulation of the tar and flgM promoters could be inferred successfully. For this analysis, we used available measurements of FliA and FlgM half-lives. This extended model contains the same number of parameters as the initial model and an improved fit between data and model is therefore not simply a consequence of increasing the degrees of freedom.

Our results also confirmed the importance of global physiological effects and the active regulation of FliA and FlgM half-lives in predicting the activity of FliA-dependent promoters. When global physiological effects were ignored, or the FliA and FlgM half-lives were set to typical values of stable E. coli proteins, a sharp drop in the quality of the prediction of the gene regulation models was observed.

CONCLUSION

More generally, under which conditions does the inclusion of the above factors lead to better results and when can they be ignored? We performed a simulation study in which we systematically varied the relative contribution of global physiological effects to cross-condition variations in the expression of a target gene and the half-lives of the regulators. These results showed that increasing half-lives of the activating transcription factor and stronger variations of global physiological effects make it more difficult to obtain good fits when using promoter activities and data on specific regulators only, respectively. While these conclusions are not surprising, it is important to emphasize that in the system studied here, where FliA and FlgM have half-lives that are exceptionally short for bacterial proteins, a considerable improvement of the fit could be obtained. For regulatory proteins with more typical half-lives, the gain may therefore be even more important than observed here.

Perspectives

Our method to more fully exploit the information contained in time-series (populationaveraged) data of the transcriptional response of bacterial cells to a changing environment depends on kinetic models of gene expression, relating the primary fluorescence and absorbance data to promoter activities and protein concentrations. In order to further improve the estimation of the biologically important quantities (promoter activity, etc.) from the primary data, we could take into account delays that are due to the maturation of GFP and the time for rounds of transcription and translation to complete. This refinement was not necessary here since the GFP reporter used in our study is fast-folding and the transcription and translation delays are very short with respect to the time-scale of the experiments.

In principle, it should be possible to apply the same approach to the inference of regulatory networks from high-throughput transcriptome data, such as DNA microarray and RNAseq data. The primary data would directly yield mRNA measurements, eliminating one step of the data treatment. However, since our algorithm relies on the interpretation of the dynamics of the system, we would need a relatively high sampling density. Furthermore, since we need to estimate the error on the change in mRNA or protein concentrations, many replicate experiments would have to be performed.

Combining a high sampling density with numerous replicates may rapidly become very 4.2 Perspectives costly. Future improvements, and the associated cost reduction, of sequencing techniques may soon make RNAseq data available for being analyzed by our algorithm.

We have found that an estimate of the protein concentration is crucial for improving the reliability of the network reconstruction. Since the measurement of mRNA is much easier than the direct measurement of protein concentrations, this implies that we need to know (or estimate) the protein half-lives. As pointed out above, in fast-growing cultures, the effective protein half-life is dominated by growth dilution. Experiments carried out in these conditions are directly amenable to analysis by our method. However, often this is not possible, for example, when measuring the transition from exponential growth to stationary phase. In this case, we would need to measure protein half-lives, for example using Western blots. However, such experiments are time and money consuming. An alternative could be to use translational fusions of the genes in the network to different flavors of fluorescent proteins. The stability of the protein would be directly correlated with the easily observable fluorescence signal. Control experiments would have to ascertain that the GFP-tag does not affect the half-life of the protein. The ideal solution, a direct observation of the proteins in the cell, for example by quantitative proteomics, remains too time and money consuming even more so than in the case of high-throughput transcriptomics.

The method we propose has the advantage that it can be used to monitor the dynamics of gene expression and global physiological effects in real time, without any additional preparation steps. However, reporter constructs have to be constructed for the genes of interest on plasmids or integrated into the chromosome. Such cloning tasks become increasingly automated and it seems possible that in the near future, many hundreds of such constructs can be assembled in parallel in a liquid handling robot. Efficient DNA assembly techniques, such as the Gibson assembly, will be optimized to streamline vector construction. We therefore anticipate that the reporter gene technology will be easily extended to other bacteria and certain eucaryotes, such as yeast.

Our approach to network reconstruction is based on the analysis of time-series data to first pinpoint possible network structures and then adjust model parameters in order to obtain a quantitative fit of the model to the experimental data. In the first step, any additional information limiting the possible network topologies could potentially be incorporated in order to reduce the search space of the algorithm. This The plots correspond to the five experimental conditions considered (∆rpoS -M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB, and WT-LB) as well as the pooling of the data sets from all five conditions. The dot in the center of each region in the plots corresponds to a tested combination of half-lives of FliA and FlgM, and thus to specific protein concentration profiles computed from the kinetic model of gene expression (Section 3.2.4). The minimal regulatory patterns were obtained by applying the minimal sign pattern algorithm (Porreca et al., 2010a) 2 was fit to the data using the promoter activity for tar (Figure 3.3) and concentrations of FliA and FlgM reconstructed from the activities of their promoters for physiologically realistic half-lives (Figure 3.12 and Figure 3.13). The parameters were estimated using a multistart global optimization algorithm (see Section 3.2.5 for details). The best fit is shown, for measured half-lives of FliA and FlgM of 30 min and 18 min, respectively (solid line, Q = 36, (k 0 , k 1 , n, θ, K) = (13.6, 206.8, 1.9, 3985, 223700)).

C. Computation of active FliA

The active regulator in the FliA-FlgM module is free FliA, that is, FliA not bound to FlgM. The active concentration of FliA can be computed from the total concentration of FliA using Eq. 2 in the main text, given a value for the equilibrium constant K

E. Additional information on plasmid construction

In order to account for the global physiological effects we used the vector from the library developed at the Weizmann Institute [START_REF] Zaslaver | A comprehensive library of fluorescent transcriptional reporters for Escherichia coli[END_REF] and we constructed a reporter for the constitutive promoter pRM of the phage lambda. The pRM promoter region was cloned into the pUA66gfp plasmid backbone using the Gibson Assembly method [START_REF] Gibson | Enzymatic assembly of overlapping DNA fragments[END_REF] and the primer sequences detailed in the table below (Table E.1). Table E.1: Primers used for the construction of pUA66pRM-gfp plasmid. The pUA66pRM-gfp plasmid was constructed with the Gibson Assembly method [START_REF] Gibson | Enzymatic assembly of overlapping DNA fragments[END_REF]. The pUA66gfp plasmid backbone was amplified using the primers gfp-fw and gfp-re. The pRM promoter region was amplified from the pZE1RMgfp plasmid Berthoumieux et al. (2013b) using the primers pRM-fw and pRM-re. pRM-fw and pRM-re contain the XhoI and BamHI restriction sites, respectively, allowing the insertion of the amplified DNA between these two sites on the pUA66gfp plasmid.

Plasmid

I

  would like to thank Delphine Ropers for all scientific discussions and practical advices during my PhD. Many thanks go to Corinne PINEL, Stéphane PINHAL, Julien DEMOL, Stéphan LACOUR and Omaya DUDIN for their availability and for helping me with (not so obvious) biological experiments. I am also grateful to Françoise de CONNINCK, Valentin ZULKOWER, Edith GRAC, Nils GIORDANO, the students of the IBIS team, the MISTIS team and the NANO-D team in INRIA for the scientific exchanges and for the pleasant moments spent at work and outside work. Finally, I express my very profound gratitude to my fiancé, Bruno Roberto FRAN-CISCATTO, for providing me with unfailing support and constant encouragements throughout the years of my PhD. I also thank my heart friends, Aline, Flora, Vincent, Bertrand, Vitor and Alina for their backing and care during both good and difficult moments. My parents receive my warmest thanks as they have always been my enlightenment. Most important, they love me unconditionally. And they introduced me to science matters. iii This work was supported by the Rhône-Alpes region (cluster ISLE, PhD grant to DS), the Investissements d'Avenir Bio-informatique programme under project RESET (ANR-11-BINF-0005), the INRIA/INSERM project ColAge, and the Agence Nationale de la Recherche under project GeMCo (ANR-2010-BLAN-0201-02).

  Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure1.2: The hierarchical transcriptional network controlling flagellar assembly and the synthesis of the chemotaxis system(Karlinsey et al., 2000a). Class 1 genes encode the proteins FlhD and FlhC, transcriptional activators that are required for expression of all the other flagellar genes. Class 2 genes encode structural components of the Hook-Basal Body (HBB) complex, as well as regulatory proteins, including the sigma factor σ 28 and the anti-σ 28 factor FlgM. σ 28 is required for the expression of all class 3 gene promoters, including those for flagellin and those related to chemotaxis and motility. FlgM binds to and inhibits σ 28 until completion of the HBB. FlgM is then secreted from the cell, and σ 28 -dependent (class 3) gene expression is initiated.

  The clustering algorithm developed by[START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF] has been applied to a variety of systems such as in[START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization[END_REF] to identify cell cycle-regulated genes of the yeast Saccharomyces cerevisiae from microarray data, and in[START_REF] Bansal | How to infer gene networks from expression profiles[END_REF] to recover network structure from steady-state and time-series data by both in silico and experimental data analysis. Although the accuracy in identifying correct interactions was fairly low, based on a large dataset (S. cerevisiae steady-state dataset, seeBansal et al.

  . Reverse engineering a Boolean network means finding a Boolean function for each gene in the network such that the observed discrete data are explained by the logics of the model. REVEAL (REVerse Engineering ALgorithm, Liang et al. (1998)) is one of the various algorithms used for the inference of Boolean networks. REVEAL uses mutual information to identify a reduced set of inputs that describe the activity of an output gene and then determines the Boolean interaction functions from the data. The algorithm was shown to perform well on a network of 50 genes, each gene having up to 3 regulators, when using for the analysis only the reduced state transitions collection (100 out of 10 15 ). At a closer look, Boolean networks cannot fully describe real gene expression profiles by their binary abstraction, therefore Boolean networks are inherently limited. As well, there are important behaviors that cannot be modeled using Boolean framework, such as amplification or addition of signals. Probabilistic Boolean networks (Akutsu 2.2 Inference of gene regulatory networks

  de[START_REF] Jong | Modeling and simulation of genetic regulatory systems: a literature review[END_REF];Porreca et al. (2010a);[START_REF] Szederkényi | Inference of complex biological networks: distinguishability issues and optimization-based solutions[END_REF];[START_REF] Yang | An analytical rate expression for the kinetics of gene transcription mediated by dimeric transcription factors[END_REF]).Second, it requires the estimation of θ i for each gene i. If the functional form of the model is not fixed, structure and parameter estimation are intertwined problems, and different methods try to isolate the two with different expedients, see Section [subsection on sign pattern] for a specific example.

  the best known algorithms for ODE network model identification are Network Identification by multiple Regression (NIR), Microarray Network Identification (MNI) and Time-Series Network Identification (TSNI) Bansal et al. (2006); Cantone et al. (2009); di Bernardo et al. (2005); Gardner et al. (2003). NIR and MNI analyse steadystate mRNA measurements, whereas TSNI uses time-series datasets.

σ

  70 -dependent transcription of the class 2 operons. The latter encode the proteins making up the flagellar motor structure as well as a major regulator of the class 3 operons, the sigma factor FliA (σ 28 ). When bound to core RNA polymerase, FliA directs the transcription of a total of 5 class 3 operons[START_REF] Keseler | EcoCyc: fusing model organism databases with systems biology[END_REF], which code for the proteins forming the filament structure of the flagellum and the chemotaxis sensing system. The aspartate chemoreceptor Tar is an example of such a class 3 protein. The action of FliA is counteracted by the anti-sigma factor FlgM, which binds to FliA and thus prevents its association with RNA polymerase. FlgM is encoded by the gene flgM, which is transcribed from both a class 2 promoter and a class 3 promoter. FlgM can be excreted from the cell through the center of the basal-body structure of the flagellum (Figure3.1).

Figure 3

 3 Figure 3.1: FliA-FlgM module. A:The regulatory circuit composed of the flagellarspecific transcription factor FliA, a sigma factor also known as σ 28 , and the anti-sigma factor FlgM forms a check-point in the transcriptional hierarchy of the motility genes in E. coli. While fliA is transcribed from a single class 2 promoter (pfliA), flgM is transcribed from both a class 2 and a class 3 promoter (pflgA and pflgM, respectively). FliA binds to RNA polymerase core enzyme and directs transcription from a total of 5 class 3 promoters[START_REF] Keseler | EcoCyc: fusing model organism databases with systems biology[END_REF], including ptar and pflgM. When bound to FlgM, FliA cannot activate transcription. When the hook basal-body (HBB) structure is in place, however, FlgM is exported from the cell, thus releasing FliA from the inactive complex. FliA is subject to proteolysis by Lon, but FlgM-binding protects FliA from degradation. The fliA promoter is auto-regulated by FliA and by a number of other regulators, most importantly the motility master regulator FlhDC. The expression of FlhDC itself is under the control of a variety of regulatory factors, including RpoS, CpxR, CsgD and Crp • cAMP. The activity of the genes is measured by fusion of their promoters to a gfp reporter gene on a low-copy plasmid. Genes are shown in grey or green and their promoter regions in red. Regulatory interactions are represented by dashed lines, association and dissociation of FliA and FlgM as well as degradation and export by solid lines. The figure does not explicitly show that fliA, flgM, and tar are included in larger transcriptional units, the fliAZY, flgAMN, flgMN and tar-tap-cheRBYZ operons[START_REF] Keseler | EcoCyc: fusing model organism databases with systems biology[END_REF]. B: Pattern of known regulatory interactions for the class 3 genes tar and flgM.

Figure 3 . 2 :

 32 Figure 3.2: Primary data and promoter activities. A: Absorbance (•, black) and fluorescence (•, blue) data, corrected for background intensities, obtained with the ∆cpxR strain transformed with the ptar -gfp reporter plasmid, grown in M9 with glucose. B: Activity of the tar promoter (•, blue), computed from the primary data as described in Section 3.2.2 and in Chapter 2. The solid line corresponds to the mean of 6 replicate absorbance measurements and the shaded region to the mean of the promoter activities ± twice the standard error of the mean.

Figure 3 . 3 :

 33 Figure 3.3: Promoter activities of genes in the FliA-FlgM module. The promoter activities of fliA (green), flgM (red), and tar (blue) measured in all experimental conditions considered in this study: ∆rpoS strain grown in M9 (∆rpoS -M9), ∆cpxR strain grown in M9 (∆cpxR-M9), ∆csgD strain grown in M9 (∆csgD-M9), ∆csgD strain grown in LB (∆csgD-LB), and wild-type strain grown in LB (WT-LB). The promoter activities have been derived from the primary data as illustrated in Figure3.2.

  2.4. Notice that the model is in agreement with the expected pattern of regulatory interactions (Figure 3.1B ). How well does this model fit the data when the total concentrations of FliA and FlgM in Eq. 3.2, p A and p M , are replaced by the measured activities of fliA and flgM, respectively? We estimated the values of the kinetic parameters c = (k 0 , k 1 , n, θ, K)

3 .

 3 INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA value of the threshold θ is similar with values of the fliA activity over all conditions.

Figure 3 . 5 :

 35 Figure 3.5: Fits of regulation function of tar to reporter gene data when replacing protein concentrations by promoter activities. The regulation function of Equations 3.1-3.2 was fit using the promoter activities for tar, fliA, and flgM shown in Figure3.3, where the latter two replace the concentrations of FliA and FlgM, respectively. The parameters were estimated using a multistart global optimization algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 33.6 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) =(7.6, 853, 1, 662, 14615). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 . 6 :

 36 Figure 3.6: Fits of regulation function of flgM to reporter gene data when replacing protein concentrations by promoter activities. The regulation function of Equations 3.1-3.2 was fit using the promoter activities for fliA, and flgM shown in Figure3.3, where the latter two replace the concentrations of FliA and FlgM, respectively. The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 23 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) =(9, 582, 1, 221, 7307). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 . 7 :

 37 Figure 3.7: FliA-FlgM module extended with activity of the gene expression machinery. A: The network is the same as in Figure 3.1, but the regulation of the motility genes by global physiological effects, in particular the activity of the gene expression machinery, has been included. These regulatory interactions are shown by bold, dashed lines. B: Pattern of regulatory interactions for the class 3 genes tar and flgM.

Figure 3 . 8 :

 38 Figure 3.8: Activities of constitutive phage promoter. The activities of the phage λ promoter pRM (black) and the activity of tar (blue) measured in all experimental conditions considered in this study. The tar promoter activities are the same as shown in Figure 3.3.

Figure 3 . 9 :Figure 3 .

 393 Figure 3.9: Minimal sign patterns for the regulation of tar and flgM when replacing protein concentrations by promoter activities and including global physiological effects. For the regulation of tar gene, the expected sign pattern (fliA, flgM, tar, pRM) = (1,-1, 0, 1) is found to be consistent with the data, e.g. the promoter is activated by the gene expression machinery and FliA and repressed by FlgM. For the regulation of flgM gene, the expected sign pattern (1,-1, 0, 1) is found to be inconsistent with the data. Similarly to results presented in Figure3.4, black arcs represent the minimal consistent sign patterns. Every consistent pattern can be obtained from one of these minimal sign patterns by turning some gray arcs into black arcs with either a line (inhibition) or an arrow (activation) end.

Figure 3 .

 3 Figure 3.11: Fits of regulation function of flgM to reporter gene data when replacing protein concentrations by promoter activities and including global physiological effects. The regulation function of Equations 3.2 -3.3 was fit using the promoter activities for fliA, and flgM shown inFigure 3.3, where the latter two replace the concentrations of FliA and FlgM, respectively. Moreover, global physiological effects are quantified by the activity of the constitutively expressed pRM promoter (Figure3.8). The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). The best fit (thick solid blue line) returns the value Q = 20 for the objective function, for the parameter vector (k 0 , k 1 , n, θ, K) =(18, 14, 1.2, 817, 7307). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 .

 3 Figure 3.12 and Figure 3.13. As can be seen, the transcriptional pulse in exponential phase (Figure3.3), leading to a transient accumulation of mRNA, is seen to be followed by the prolonged presence of stable protein, indicating that the promoter activity may indeed not be a good proxy for the protein concentration. Unfortunately, reporter concentrations are not always representative of the concentrations of proteins of interest, that is, proteins naturally expressed from a promoter. Post-transcriptional regulation and coding bias may cause divergent synthesis rates. The main bias, however, comes from the fact that the two proteins generally have different half-lives and thus different degradation rates (de[START_REF] Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF].Available data in the literature indicate that the half-lives of FliA and FlgM are much shorter than the 19 h of the GFP reporter. The measured half-lives of FliA and

Figure 3 .Figure 3 .

 33 Figure 3.13: Estimates of FlgM concentrations from reporter gene data. Concentrations of FlgM (dashed line) computed from the flgM promoter activity (solid line) in all experimental conditions considered in this study. The flgM activities are the same as shown in Figure3.3. The dark red line represents the concentration of the reporter protein, while the light red line represents the reconstructed concentration for the measured half-live of 18 min. Promoter activity has been normalized with respect to the maximum of the upper limit of its confidence interval in each condition. All protein concentrations have been normalized with respect to the maximum of the upper limit of the confidence interval of the reporter concentration in each condition. The shaded region corresponds to the mean of the promoter activities ± twice the standard error of the mean.

Figure 3 .

 3 Figure 3.15: Minimal patterns of regulatory interactions for flgM over a range of physiologically realistic half-lives. The minimal regulatory patterns for the gene flgM in the motility network of Figure 3.7 as a function of the half-lives of FliA and FlgM.Similarly to Figure3.14, the plots correspond to the four of the experimental conditions considered (∆rpoS -M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB) as well as the pooling of the data sets from all five conditions. The condition WT-LB was not used in the analysis of the regulation of the flgM promoter. The dot in the center of each region in the plots corresponds to a tested combination of half-lives of FliA and FlgM, and thus to specific protein concentration profiles computed from the kinetic model of gene expression (Section 3.2.4). The minimal regulatory patterns were obtained by applying the minimal sign pattern algorithm(Porreca et al., 2010a). The color codes represent the different categories of minimal signal patterns inferred. A region is colored green if the expected regulatory patterns is among the minimal sign patterns returned by the algorithm, and yellow if it is compatible with the returned sign patterns. A region is colored red if none of the returned sign patterns is consistent with the data. Two examples of inconsistent sign patterns are shown. The values of the half-lives are represented in log.

Figure 3 .

 3 Figure 3.16: Fits of regulation function of tar to reporter gene data when reconstructing protein concentrations from the reporter gene data and including global physiological effects. The regulation function of Eqs. 3.2-3.3 was fit to the data using the promoter activity for tar (Figure3.3), concentrations of FliA and FlgM reconstructed from the activities of their promoters for physiologically realistic half-lives (Figure3.12 and Figure3.13), and the activity of the constitutively expressed pRM promoter quantifying global physiological effects (Figure3.8). The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). Three fits are shown, namely the best fit for measured half-lives of FliA and FlgM of 30 min and 18 min, respectively (solid line, Q = 24.7, (k 0 , k 1 , n, θ, K) = (0.3, 4.6, 2.4, 3030, 223750)) and two other fits for comparable half-lives (dashed lines, Q = 24.1, (k 0 , k 1 , n, θ, K) = (0.2, 4.7, 2.2, 4535, 222800) and Q = 25, (k 0 , k 1 , n, θ, K) = (0.3, 4.6, 2.4, 2800, 162000) ). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 .

 3 Figure 3.17: Fits of regulation function of flgM to reporter gene data when reconstructing protein concentrations from the reporter gene data and including global physiological effects. The regulation function of Eqs. 3.2-3.3 was fit to the data using the promoter activity for flgM (Figure3.3), concentrations of FliA and FlgM reconstructed from the activities of their promoters for physiologically realistic half-lives (Figure3.12 and Figure3.13), and the activity of the constitutively expressed pRM promoter quantifying global physiological effects (Figure3.8). The parameters were estimated using a hybrid genetic algorithm (see Section 3.2.5 for details). Three fits are shown, namely the best fit for measured half-lives of FliA and FlgM of 30 min and 18 min, respectively (solid line, Q = 27.7, (k 0 , k 1 , n, θ, K) = (0.3, 6.1, 2, 3170, 223750)) and two other fits for comparable half-lives (dashed lines, Q = 26.3, (k 0 , k 1 , n, θ, K) = (0.4, 5.9, 2.3, 2358, 279500) andQ = 27.1, (k 0 , k 1 , n, θ, K) = (0.3, 6, 2.1, 2760, 162000) ). The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 .

 3 Figure 3.14 shows the results for the structural inference of tar regulators. As can be seen, almost all combinations of half-lives are compatible with activation of tar by FliA and the gene expression machinery and with inhibition by FlgM. This means that the structure of interactions is robust over a range of half-lives, a desirable property for network inference.Figure 3.18 illustrates that the quantitative regulation function of tar activity obtained is more precise than in all other previously considered situations (Q = 20.9), while the parameter values are similar to those obtained in the previous sections. Although we substantially relaxed the possible half-live values of FliA and FlgM, it is remarkable that the optimal values are close to the reported values for LB medium (Figure 3.18). This emphasizes the importance of active degradation of FliA and secretion of FlgM for the dynamics of the motility network. Moreover, while the proportion of FliA released by FlgM varies across conditions, most FliA is predicted to be free over the duration of the experiment (Figures C.1, C.2, C.3, C.4 in Appendix C).

  Figure D.1-F in Appendix D). The 3.1 Results half-lives obtained in the case of the best fit vary around the measured values of halflives for FliA and FlgM. We obtain a more stable half life for FliA (50 min) and FlgM (27 min) in the ∆rpoS -M9 condition. The values of half-live in all conditions can be seen in the legend of Figure 3.18.

Figure 3 .

 3 Figure 3.19: Fits of regulation function of flgM to reporter gene data when reconstructing protein concentrations from the reporter gene data for physiologically realistic half-lives and including global physiological effects. Similarly to Figure3.18 the half-lives have been estimated from the data, within a physiologically plausible range. Two examples of fits are shown, namely the best fit for estimated half-lives of FliA and FlgM (solid line, Q = 25, (k 0 , k 1 , n, θ, K) = (0.45, 5.9, 2.4, 2930, 222850)) and another example of a high-ranking fit (dashed line, Q = 25.4, (k 0 , k 1 , n, θ, K) = (0.4, 6.1, 2.3, 3030, 279550)). In the case of the best fit obtained, the half-lives of FliA are equal to (24 min, 35 min, 27 min, 50 min) in the (∆rpoS, ∆cpxR, ∆csgD-M9, ∆csgD-LB) conditions, respectively. The half-lives of FlgM are equal to (27 min, 11 min, 11 min, 9 min) in the (∆rpoS, ∆cpxR, ∆csgD-M9, ∆csgD-LB) conditions, respectively. In the case of the second fit, the half-lives of FliA are equal to (24 min, 35 min, 27 min, 27 h) and the half-lives of FlgM are equal to (18 min, 13 min, 20 min, 27 min) in the above precised experimental conditions, respectively. The mean of the promoter activity of tar (thin solid blue line) and confidence intervals (shaded blue regions) are also shown in the figure.

Figure 3 .

 3 Figure 3.20: Heatmap of the fitting residuals, given by the value of the objective function Q, for simulated data generated for different protein half-lives and for different strengths of global physiological effects. A: For all different combinations of 33 half-lives of FlgM (horizontal axis) and FliA (vertical axis), the residual of the fit by a model ignoring protein kinetics is represented by the color code reported in the right bar.The combination corresponding to the measured half-lives in rich LB medium is marked with a light blue square (18 min for FlgM, 30 min for FliA). B: For 26 different values of the strength parameter α, defined in Eq. 3.4, the residual of the fit by a model ignoring global physiological effects is represented by the color code. The value corresponding to the real data is marked with a light blue rectangle (α = 1).

Figure 3 .

 3 21 gives an example of background correction of absorbance and fluorescence data, in the case of the tar reporter in the ∆cpxR mutant strain.

Figure 3 .

 3 Figure 3.21: Illustration of data analysis procedures. Absorbance and fluorescence data acquired with the ∆cpxR mutant strain carying a pUA66tar-gfp plasmid, grown in M9 with glucose. A: Primary (uncorrected) absorbance (•, grey), background absorbance (•, red), and corrected absorbance (•, black). B: Calibration curve obtained by measuring the autofluorescence of the wild-type strain without plasmid. Primary fluorescence data are plotted against (corrected) absorbance data and the curve is obtained by fitting a smoothing spline. C: Primary fluorescence data (•, grey), and the corrected fluorescence (•, blue) obtained after subtracting the fluorescence of the background (•, red) as in Eq. 3.6. D: Promoter activity of tar (•, blue) computed from the corrected absorbance (-, grey) and corrected fluorescence by means of Eq. 3.10. E: Protein concentration of tar (•, blue) computed for a half-life of 2 h from the corrected absorbance (-, grey) and corrected fluorescence measurements using Eq. 3.7. F: Protein concentration of tar (•, blue) computed for a half-life of 18 h from the corrected absorbance (-, grey) and corrected fluorescence measurements using Eq. 3.7.

Figure 3 .

 3 22 illustrates the effects of variable half-lives and extrapolation of promoter activities on the computation of FliA and FlgM concentrations in a WT strain. Moreover, in various experimental conditions (rich medium) the activity of the promoters can only be observed when the fluorescence 3.2 Methods and materials intensity overreaches the value of background fluorescence. When fluorescence intensity of FliA and FlgM was not observable before actual expression of genes we have assumed it to be 0 and we have interpolated the promoter activity values towards 0 for this part of the experiment.Figure 3.22 shows an example of effects of variable half-lives and interpolation of promoter activities on the computation of protein concentrations of FlgM and FliA in a WT strain.

  2. A set of measurementsD = {(x k , f k ) : k = 1, . . . ,m} of the concentration vectors x and the corresponding target promoter activities f are available, along with confidence intervals f k ± k and xk ± e k .

Figure 3 .

 3 Figure 3.22: Effect of variable half-lives and promoter activity extrapolation on initial protein concentrations in a WT strain. A: The observed promoter activtity of flgM (dashed line, red) and its etrapolation (*, red). B : The observed promoter activtity of fliA (dashed line, green) and its etrapolation (*, green). C : The effect on the computation of the protein concentration of FlgM when taking into account the extrapolation of its promoter activity in A and the initial condition. Protein concentration was computed for an initial half-life of 3 h and a short half life of 18 min. D: The effect on the computation of the protein concentration of FliA when taking into account the extrapolation of its promoter activity in B and the initial condition. Protein concentration was computed for an initial half-life of 2 h and a short half life of 30 min. The hatched regions correspond to the regions where the activity of promoters was extrapolated, i.e. [0, 410 min] for flgM and [0, 340 min] for fliA. The promoter activities have been derived from the primary data as illustrated in Figure3.2. The shaded regions correspond to the mean of the promoter activities and protein concentrations, respectively, ± twice the standard error of the mean. The absorbance is drawn in solid, grey lines.

  and π j = 0 if f is independent of x j , j = 1, . . . , n. The sign pattern encodes the directed, signed graph of the regulation of the gene under consideration by all possible regulators in the net-3.2 Methods and materialswork (compare Figure3.23 B-E). As for assumption 2, data may come from several gene reporter experimental scenarios (different strains and media) and is provided in the required form by the processing of the previous Section 3.2.2, where (x k , f k ) is the measurement average at time t k , while e k and k are fixed to twice the standard error of the mean (x k , f k ). Also observe that dependence of confidence intervals on index k is explicitly taken into account.

For

  any data point (x, f ) ∈ D, let (x, f ) and (x, f ) indicate the confidence bounds f = f + and f = f -, in the same order, and similarly xj = xj + e and xj = xj -e, with j = 1, . . . , n. Let the complexity c of a sign pattern π be the number of nonzero 3. INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA entries of π. The algorithm is divided into two phases, conceptually organized as follows (see Figure 3.23 for reference). Computation of inconsistent patterns Π from data D • Set Π = ∅

  14, and is analogous for flgM. The regressors take different forms in consecutive sections of this thesis, consisting either of the activities f A and f M of the fliA and flgM promoters (x = (f A , f M )) or the reconstructed concentrations p A and p M of FliA and FlgM (x = (p A , p M )). In all sections, c = (k 0 , k 1 , n, θ, K). The superscript symbol s indicates the experimental condition, where s ∈ S = {∆rpoS-M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB, WT-LB}.

  transcriptional and translational machinery of the cell, as well as other intrinsic global physiological effects, may induce major changes in gene expression over the course of an experiment . As[START_REF] Lovén | Revisiting global gene expression analysis[END_REF] point out, many transcriptome studies assume that the total quantity of RNA is similar between different experimental conditions and use this quantity for normalization of the data. As a consequence, a global increase or decrease of transcriptional activity across conditions may lead to erroneous interpretation of the experimental data and the inference of spurious regulatory interactions. In this thesis, we have developed an improved model that, in addition to specific regulatory interactions, accounts for global regulatory effects of bacterial cells. We estimate the global physiological state of the cell from the activity of a constitutively expressed gene and whose expression only depends on the activity of the transcriptional and translational machinery.Our analysis of reporter gene datasets thus makes it possible to account for the difference between mRNA and protein concentrations as well as for global physiological effects. We have validated the adequacy of the models to describe bacterial gene expression by comparing the deduced promoter activities with independent RT-qPCR measurements. We have shown that the inclusion of information about both global physiological effects and protein concentrations can improve the inference of regulatory interactions and the identification of quantitative regulation functions from time-series data. Compared to the classical inference approach, i.e., the inference of structure and 4.1 Summary of results quantitative functions of gene regulatory networks from promoter activities, the inclusions of global regulatory effects significantly improves the prediction results. A further improvement is achieved by explicitly considering proteins as the regulators by taking into their half-lives.

A.

  Figure A.1: The flgAMN operon.

FigureFigure A. 2 :

 2 Figure A.2.

Figure B. 3 :

 3 Figure B.2: Minimal patterns of regulatory interactions for flgM over a range of physiologically realistic half-lives. The minimal regulatory patterns for the gene flgM in the motility network of Figure 3.7 as a function of the half-lives of FliA and FlgM.Similarly to Figure3.14, the plots correspond to the four of the experimental conditions considered (∆rpoS -M9, ∆cpxR-M9, ∆csgD-M9, ∆csgD-LB) as well as the pooling of the data sets from all five conditions. The condition WT-LB was not used in the analysis of the regulation of the flgM promoter. The dot in the center of each region in the plots corresponds to a tested combination of half-lives of FliA and FlgM, and thus to specific protein concentration profiles computed from the kinetic model of gene expression (Section 3.2.4). The minimal regulatory patterns were obtained by applying the minimal sign pattern algorithm(Porreca et al., 2010a). The color codes represent the different categories of minimal signal patterns inferred. A region is colored green if the expected regulatory patterns is among the minimal sign patterns returned by the algorithm, and yellow if it is compatible with the returned sign patterns. A region is colored red if none of the returned sign patterns is consistent with the data. Two examples of inconsistent sign patterns are shown. The values of the half-lives are represented in log.

Figure B. 4 :Figure B. 5 :Figure B. 6 :

 456 Figure B.4: Fits of regulation function of flgM to reporter gene data when reconstructing protein concentrations from the reporter gene data and ignoring global physiological effects. The regulation function of Eqs. 3.1-3.2 was fit to the data using the promoter activity for flgM (Figure3.3) and concentrations of FliA and FlgM reconstructed from the activities of their promoters for physiologically realistic halflives (Figure3.12 and Figure3.13). The parameters were estimated using a multistart global optimization algorithm (see Section 3.2.5 for details). The best fit is shown, for measured half-lives of FliA and FlgM of 30 min and 18 min, respectively (solid line, Q = 27, (k 0 , k 1 , n, θ, K) =(31.2, 246.8, 3, 2353, 223700)).

FigureFigureFigure C. 3 :

 3 Figure C.1: Free and total concentration of FliA when using promoter activities. The concentration of free FliA (solid line, green) is computed by means of Eq. 3.2 in Chapter 3, for the optimal fit shown in Figure3.5 in Chapter 3. The shaded regions represent the confidence intervals of total FliA and correspond to the mean of the promoter activities for 6 replicates ± twice the standard error of the mean. The threshold parameter θ is shown as a dashed green line.

  pRM-fw: GAGGC CCTTT CGTCT TCACC TCGAG CCTAT CACCG CCAGA pRM-re: TTCTT AAATC TAGAG GATCC GGTTT CTTTT TTGTG CTGAT gfp-fw: ATCAG CACAA AAAAG AAACC GGATC CTCTA GATTT AAGAA gfp-re: TCTGG CGGTG ATAGG CTCGA GGTGA AGACG AAAGG GCCTC

  CONTENTS 2.3 Reporter gene data analysis . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Measurement models . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Constitutive promoters . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Reconstruction of promoter activities . . . . . . . . . . . . . . . . 2.3.5 Reconstruction of protein concentrations . . . . . . . . . . . . . . Derivation of regulation function of motility genes . . . . . . . . 3.2.5 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . 3.2.6 Validation of reporter gene data using quantitative RT-PCR . .

	3 Inference of quantitative models of bacterial promoters from time-
	series gene expression data

3.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Monitoring the transcriptional response of the FliA-FlgM module 3.1.2 Identification of gene regulation functions from promoter activities 3.1.3 Identification of gene regulation functions from promoter activities including global physiological effects . . . . . . . . . . . . . . 3.1.4 Identification of gene regulation functions from estimates of protein concentrations . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.5 Determination of conditions in which protein half-lives and global physiological effects are important . . . . . . . . . . . . . . . . . 3.2 Methods and materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Strains and growth conditions . . . . . . . . . . . . . . . . . . . . 3.2.2 Experimental monitoring of gene expression in real time and data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2.1 Background subtraction . . . . . . . . . . . . . . . . . . 3.2.2.2 Computation of promoter activity and protein concentrations . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Computation of minimal consistent sign patterns of regulatory interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2.4

  . However, identification of nonlinear ODE parameters often results in non-convex problems, revealing a series of difficulties if standard local optimization methods are used, such as converging to local solutions or bad scalability for large systems. Global optimization methods have been devel-

oped to seek globally optimal solutions (see

[START_REF] Floudas | A review of recent advances in global optimization[END_REF]

for a review;

[START_REF] Banga | Optimization in computational systems biology[END_REF]

;

[START_REF] Banga | Dynamic optimization of bioprocesses: efficient and robust numerical strategies[END_REF]

;

Vilas et al. (

  STATE OF THE ARTand by replacing k p by k p (t), and k r by k r (t). Analogously to Eq. 2.19, the generalized expression for the synthesis rate of the reporter protein becomes:f (t) = k m kr (t) g global (t) g specif ic (t), (2.22)which is decomposed in a part due to the activity of the gene expression machinery (k m kp (t) g global (t)) and a part due to specific effects of transcription regulators (g specif ic (t)). By the same reasoning as for the classic measurement models, this expression remains proportional to the synthesis rate of the protein of interest (with

by positing

g(t) = k m g global (t) g specif ic (t),

(2.21)

2.

  .1).

	Strain	Characteristics	Reference or source
	WT	E. coli BW25113	Baba et al. (2006)
	WTpRM E. coli BW25113 pRM-gfp::intS	This study
	∆rpoS	E. coli BW25113 ∆rpoS	Dudin et al. (2013)
	∆cpxR	E. coli BW25113 ∆cpxR	Dudin et al. (2013)
	∆csgD	E. coli BW25113 ∆csgD	Dudin et al. (2013)

Table 3

 3

.1: Strains used in this study.

Table 3 .

 3 2: Plasmids used in this study.

	3. INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL
	PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA
	Plasmid	Characteristics	Reference or source
	pUA66gfp	Kan r , pSC101ori, gfpmut2	Zaslaver et al. (2006)
	pUA66fliA-gfp	Kan r , pSC101ori , fliA -gfpmut2	Zaslaver et al. (2006)
	pUA66flgM-gfp Kan r , pSC101ori , flgM -gfpmut2 Zaslaver et al. (2006)
	pUA66flgA-gfp	Kan r , pSC101ori , flgA -gfpmut2 Zaslaver et al. (2006)
	pUA66tar-gfp	Kan r , pSC101ori , tar -gfpmut2	Zaslaver et al. (2006)
	pUA66pRM-gfp Kan r , pSC101ori , pRM -gfpmut2	This study

  are the measured C T values for gfp and tar, respectively, ∆C gf p T

		.16)
	where C gf p T	and C tar T

  . The color codes represent the different categories of minimal signal patterns inferred. A region is colored green if the expected regulatory patterns is among the minimal sign patterns returned by the algorithm, and yellow if it is compatible with the returned sign patterns. A region is colored red if none of the returned sign patterns is consistent with the data. Two examples of inconsistent sign patterns are shown. The values of the half-lives are represented in log.

With the definition of the initial protein concentration and additional information on the half-life of the protein (degradation constant γ p ) p(t) can be computed by numerical solution of the above ODE.

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

relevant parameters in the model (k 0 , k 1 , n, θ, K) were fixed in agreement with the best fit obtained for the reference half-lives of 30 min for FliA and 18 min for FlgM, shown in Figure 3.16. More precisely, the relative position of the parameter values within the interval of physiologically plausible values, which may depend on the FliA and FlgM concentrations, as explained in Section 3.2, was conserved across conditions. Activity profiles of tar were then generated in accordance with Eqs. 3.2-3.3 based on the experimentally measured pRM activities. We then attempted to identify from these simulated data a gene regulation model accounting for the global physiological effects, but using promoter activities in place of FliA and FlgM concentrations. The results are reported in Figure 3.20.

As can be seen, the quality of the fit decreases with longer half-lives of FliA, but is rather insensitive to the half-life of FlgM. The strong dependency on the half-life of FliA shows that, in general, accounting for slow protein kinetics is important, but that promoter activities can be safely used in place of protein concentrations for very fast-degrading proteins. This is intuitively explained by the fact that fast-degrading protein concentration profiles reproduce promoter activity profiles quite closely, while this is not true in case of slow degradation (Figure 3.12). The relative insensitivity to FlgM half-lives can be explained by the fact that, in the time window considered in our experimental set-up, a good fit requires most FliA to be free (Appendix C).

Longer half-lives, and therefore higher concentrations of FlgM, favor lower free FliA concentrations, but this tendency is compensated in the parameter optimization process by higher values for the equilibrium constant K. The actually measured reference half-lives of 18 min for FlgM and 30 min for FliA are located in the upper left corner of Figure 3.20A, where fitting residuals are comparably small. Therefore, for networks involving regulators with longer half-lives than the exceptionally short half-lives observed for FliA and FlgM, it will be even more critical to account for protein kinetics.

To evaluate the importance of global physiological effects, starting from the experimentally measured pRM activity profiles, we simulated global physiological effects of different strength. In particular, we rescaled the variations of f const (t) around its temporal mean across all conditions, fconst , by a factor α ranging from 0 (no variability, no regulatory effect) to 1 (measured variability, moderate regulatory effect) and 1.25 (increased variability, strong regulatory effect). That is, synthetic activity profiles of FliA-dependent promoters were generated in accordance with the model

INFERENCE OF QUANTITATIVE MODELS OF BACTERIAL PROMOTERS FROM TIME-SERIES GENE EXPRESSION DATA

In practice, the above operations can be made computationally efficient. Notably, in our implementation, "for" loops and "if" tests are replaced by suitable algebraic and Boolean matrix operations in Matlab.

Derivation of regulation function of motility genes

We develop a kinetic model for the regulation of the expression of tar as a function of the total concentrations of FliA and FlgM (see Figure 3.1 in Section 3.1.1). The model is based on a quasi-equilibrium approximation of the mass-action kinetics for the formation of the FliA•FlgM complex, and a phenomenological Hill-type regulatory law of tar expression by free FliA.

Let p A,f ree , p M,f ree and p AM denote the concentrations of free FliA, free FlgM and FliA•FlgM, respectively, and let p A , p M denote total concentrations for FliA and FlgM.

Assuming complex formation and dissociation are fast events relative to gene expression and protein degradation, we make the approximation

with k -> 0 and k + > 0. Using the facts that p A = p A,f ree + p AM and p M = p M,f ree + p AM , substitution into the above to eliminate p M,f ree and p AM from the equation yields

which is a second-order polynomial equation in p A,f ree . The solution of the equation (3.12) with K = k -/k + , which is a function of the total concentrations p A and p M [START_REF] Buchler | Molecular titration and ultrasensitivity in regulatory networks[END_REF]. Only the free FliA molecules regulate the expression of the tar promoter, and we quantify the regulatory effect by the law The figure shows the promoter activity of tar with respect to global physiological effects (•, blue). The promoter activities were derived from corrected absorbance and fluorescence, measured using plasmids expressing the GFP reporter for the phage promoter pRM and the tar promoter. The shaded regions represent confidence intervals computed as ± the standard error of the mean of 5 replicates. B :The figure reports the expression of tar gene with respect to the mRNA quantity of pRM (WTpRM strain, Section 3.2.1) measured by qPCR (•, black). We have normalized with respect to the observed Tar mRNA quantity in exponential growth following Eq. 3.16. Expression of tar gene is observed to be maximal around 700 min. The results obtained by using the two independent techniques are in good agreement. The errorbars were computed from the standard error of the mean of 4 replicates. The absorbances have been also plotted on the figures (solid line, grey).

additional information can come from any of a number of sources: biological data that exclude a particular network structure, "classical" inference methods using steadystate data of the same system to limit the possible network topologies, or any other modeling approaches or measurements that fix certain parameters of the quantitative model. We therefore consider our algorithm a further improvement to the general problem of network inference, adding a new powerful tool that can be combined with existing methods to reliably deduce the underlying regulatory structure from time-series expression data. .4: Free and total concentration of FliA when using reconstructed protein concentrations for physiologically-realistic half-lives estimated from the data, and including global physiological effects for tar regulation function. The concentration of free FliA (solid line, green) is computed by means of Eq. 3.2 in Chapter 3, for the optimal fit shown in Figure 3.18 in Chapter 3. The shaded regions represent the confidence intervals of total FliA and correspond to the mean of the promoter activities for 6 replicates ± twice the standard error of the mean. The threshold parameter θ is shown as a dashed green line.

E. ADDITIONAL INFORMATION ON PLASMID CONSTRUCTION