<figure><img src="image_1.png" /></figure>

HAL
open science

# Modelling biochemical reaction networks in bacteria – From data to models and back

Delphine Ropers

► To cite this version:

Delphine Ropers. Modelling biochemical reaction networks in bacteria – From data to models and back. Bioinformatics [q-bio.QM]. Université Claude Bernard Lyon I, 2021. tel-03252736

HAL Id: tel-03252736

https://inria.hal.science/tel-03252736

Submitted on 7 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
HABILITATION A DIRIGER DES RECHERCHES

Délivrée par :

l'Université Claude Bernard Lyon 1

Ecole Doctorale 341

E2M2 : Evolution Ecosystèmes Microbiologie Modélisation

Spécialité : Bioinformatique

Soutenue publiement le 17/05/2021, par :

Delphine Ropers

Modelling biochemical reaction networks in bacteria

From data to models and back

Devant le jury composé de :

<table><tr><td><strong>Alexander BOCKMAYR</strong><br>Professeur, Freie Universität Berlin</td><td>Rapporteur</td></tr><tr><td><strong>Matthieu JULES</strong><br>Professeur, AgroParisTech / Chargé de recherche, INRAE</td><td>Rapporteur</td></tr><tr><td><strong>Vassily HATZIMANIKATIS</strong><br>Professeur, Ecole Polytechnique Fédérale de Lausanne</td><td>Rapporteur</td></tr><tr><td><strong>Sandrine CHARLES</strong><br>Professeure, Université de Lyon I</td><td>Présidente</td></tr><tr><td><strong>Muriel COCAIGN-BOUSQUET</strong><br>Directrice de recherche, INRAE</td><td>Examinatrice</td></tr><tr><td><strong>Johannes GEISELMANN</strong><br>Professeur, Université Grenoble - Alpes</td><td>Examineur</td></tr></table>
<figure><img src="image_2.png" /></figure>
# Université Claude Bernard – LYON 1

Administrateur provisoire de l'Université

M. Frédéric FLEURY

Président du Conseil Académique

M. Hamda BEN HADID

Vice-Président du Conseil d'Administration

M. Didier REVEL

Vice-Président du Conseil des Etudes et de la Vie Universitaire

M. Philippe CHEVALLIER

Vice-Président de la Commission de Recherche

M. Jean-François MORNEX

Directeur Général des Services

M. Pierre ROLLAND

## COMPOSANTES SANTE

Département de Formation et Centre de Recherche
en Biologie Humaine

Directrice : Mme Anne-Marie SCHOTT

Faculté d'Odontologie

Doyenne : Mme Dominique SEUX

Faculté de Médecine et Maïeutique Lyon Sud - Charles Mérieux

Doyenne : Mme Carole BURILLON

Faculté de Médecine Lyon-Est

Doyen : M. Gilles RODE

Institut des Sciences et Techniques de la Réadaptation (ISTR)

Directeur : M. Xavier PERROT

Institut des Sciences Pharmaceutiques et Biologiques (ISBP)

Directrice : Mme Christine VINCIGUERRA

## COMPOSANTES & DEPARTEMENTS DE SCIENCES & TECHNOLOGIE

Département Génie Electrique et des Procédés (GEP)

Directrice : Mme Rosaria FERRIGNO

Département Informatique

Directeur : M. Behzad SHARIAT

Département Mécanique

Directeur : M. Marc BUFFAT

Ecole Supérieure de Chimie, Physique, Electronique (CPE Lyon)

Directeur : Gérard PIGNAULT

Institut de Science Financière et d'Assurances (ISFA)

Directeur : M. Nicolas LEBOISNE

Institut National du Professorat et de l'Education

Administrateur Provisoire : M. Pierre CHAREYRON

Institut Universitaire de Technologie de Lyon 1

Directeur : M. Christophe VITON

Observatoire de Lyon

Directrice : Mme Isabelle DANIEL

Polytechnique Lyon

Directeur : Emmanuel PERRIN

UFR Biosciences

Administratrice provisoire : Mme Kathrin GIESELER

UFR des Sciences et Techniques des Activités Physiques et Sportives (STAPS)

Directeur : M. Yannick VANPOULLE

UFR Faculté des Sciences

Directeur : M. Bruno ANDRIOLETTI
# Abstract

## English abstract

With the advent of new technologies, experimental data in biology has exploded in size and complexity. It is now possible to simultaneously quantify different components of the cell at metabolic, transcriptomic, proteomic, and phenotypic levels. Connecting these different multi-scale and dynamic datasets provides an integrated view of cellular growth and informs us about the underlying molecular networks of genes, RNAs, proteins and metabolites that control the adaptation of the cell to the environment. This is the perspective offered by mathematical modelling and computer simulation, allowing the association of different microscopic and macroscopic scales. This is a difficult problem however, because of the noise and the heterogeneity of the data, and of the size and the nonlinearity of the models. As a consequence, a large number of datasets are only partially analysed and underexploited. This manuscript describes the work I have carried out to improve the utilization of experimental data to gain a better understanding of the adaptation of bacterial growth to a changing environment. This work has been carried out within the Ibis project-team (Inria, Université Grenoble Alpes) with my colleagues, especially the students that I have had the chance to supervise.

After the introductory Chapter 1, I describe in Chapter 2 the modelling of cellular networks using ordinary differential equations as well as simplification and approximation of the models depending on the nature of the available data and the questions addressed. These principles are applied in Chapter 3 to the qualitative analysis of the dynamics of gene networks in the context of the carbon starvation response in *Escherichia coli* bacteria. With the general trend of biology becoming increasingly quantitative, modelling studies require obtaining reliable gene expression and metabolomic data, the analysis of which requires the development of suitable methods described in Chapter 4. Chapter 5 examines the strong link between the activity of the cellular gene expression machinery and bacterial growth rate. This understanding is used to develop a synthetic strain of *E. coli* whose growth control makes it possible to divert the flow of precursors for growth towards the bioproduction of molecules of biotechnological interest. In Chapter 6, large-scale reconstructions of central carbon metabolism are used as platforms to interpret datasets regarding the post-transcriptional regulation of central carbon metabolism in *E. coli*. Chapter 7 is dedicated to the genome-scale analysis of mRNA decay by means of dynamic transcriptomics data. I describe in Chapter 8 ongoing and future projects towards the integrative analysis of microbial growth and resource allocation strategies. The scientific developments of these projects are expected to shape my own research activity in the coming years and that of the future project-team, under creation, that I will lead.
Long French abstract

Avec l'arrivée des nouvelles technologies, les données expérimentales en biologie ont explosé en taille et complexité. Il est désormais possible de quantifier en même temps différents composants de la cellule au niveau métabolique, transcriptomique, protéomique et de caractéristiques phénotypiques comme le taux de croissance. Relier ces différents jeux de données multi-échelles et dynamiques permet d'obtenir une vision intégrée de la croissance cellulaire, en nous renseignant sur la façon dont les réseaux moléculaires sous-jacents de gènes, ARN, protéines et métabolites contrôlent l'adaptation des cellules à leur environnement. C'est le cadre qu'offrent la modélisation mathématique et la simulation informatique, en permettant d'associer les différentes échelles microscopiques et macroscopiques. C'est cependant un problème difficile, du fait du bruit et de l'hétérogénéité des données d'une part, et de la taille et la forme non-linéaire des modèles d'autre part. La conséquence est qu'un grand nombre de jeux de données ne sont que partiellement analysés et sous-exploités.

Ce manuscrit décrit les travaux que j'ai menés pour améliorer l'utilisation de données expérimentales afin d'obtenir une meilleure compréhension de l'adaptation de la croissance bactérienne à un environnement changeant. Ces travaux ont été menés au sein de l'équipe-projet Ibis (Inria, Université Grenoble Alpes) avec mes collègues, en particulier les étudiants que j'ai eu la chance d'encadrer. Après un premier chapitre d'introduction, je décris en chapitre 2 les concepts de base de la modélisation des réseaux biochimiques. Je détaillerai en particulier les reconstructions du métabolisme cellulaire à l'échelle du génome et la modélisation cinétique des réactions enzymatiques, dont les concepts sont utilisés dans plusieurs travaux présentés dans ce manuscrit. La grande dimension et non linéarité des modèles cinétiques complique l'estimation de leurs paramètres et l'analyse de leur dynamique. Je présenterai des travaux sur des simplifications appropriées pour ces modèles selon la nature des données à disposition et les questions abordées, comme la réduction de modèles d'équations différentielles ordinaires (ODE) par séparation des échelles de temps ou l'approximation des modèles ODE par des modèles linéaires par morceaux. Du fait de leur dérivation rigoureuse, les modèles simplifiés retiennent les principales caractéristiques des modèles ODE. Ces approches seront utilisées pour les différents modèles dynamiques présentés dans ce manuscrit.

Dans le chapitre 3, je présente des travaux d'analyse de la dynamique d'un réseau de régulation génique contrôlant la réponse à la privation en carbon de la bactérie *Escherichia coli*. Lors de ces travaux, l'absence de données quantitatives dans la littérature ne permettait pas d'utiliser un modèle ODE pour décrire la dynamique du système. J'ai plutôt analysé la dynamique d'une version linéaire par mocreaux de ce modèle par une approche de modélisation et simulation qualitative. Je décrirai le principe de cette approche avec un exemple simple et son application à l'étude du réseau de la réponse au manque de source de carbone. Cette approche a permis pour la première fois de relier la croissance d'E. coli avec les principaux régulateurs transcriptionnels de la bactérie, et de comprendre les cascades de régulations mises en place lors de la réponse à une privation en glucose ou du rédémarrage de croissance sur ce sucre.

L'évolution de la biologie en une science quantitative permet d'obtenir de nombreuses
données d'expression génique et du métabolisme cellulaire. La fiabilité de ces données nécessite
le développement de méthodes d'analyse adaptées décrites dans le chapitre 4. Je décrirai des
travaux sur l'analyse de données de gènes rapporteurs et l'analyse de données de métabolomique
afin de pouvoir reconstruire des profils d'activités de promoteurs et de concentrations de pro-
téines dans le premier cas, et des vitesses d'import et secrétion de métabolites extracellulaires,
ainsi que des taux de croissance dans le second cas. Les données quantitatives utilisées dans le
reste du manuscrit ont été analysées grâce à ces approches.

Le chapitre 5 s'intéresse au lien étroit entre activité de la machinerie cellulaire d'expression génique et taux de croissance bactérien. A l'aide de modèles simples intégrant des données expérimentales de gènes rapporteurs, nous montrons le rôle clé joué par la machinerie d'expression génique dans l'adaptation globale de l'expression des gènes au cours de la croissance. Ces travaux montrent que le fonctionnement des réseaux biochimiques ne peut être déconnecté de l'état physiologique de la cellule. Cette compréhension est utilisée pour l'ingénierie d'une souche d'E. coli synthétique dont le contrôle de la croissance permet de divertir les flux de précurseurs pour la croissance vers la bioproduction de molécules d'intérêt biotechnologique.

Dans le chapitre 6, de grandes reconstructions du métabolisme et différents jeux de données (métabolomique, activités spécifiques) sont utilisées pour étudier la régulation post-transcriptionnelle du métabolisme central carboné chez *E. coli*. Ces travaux ont permis d'expliquer les conséquences physiologiques de l'atténuation du gène de la protéine CsrA et d'identifier des ARNm cibles de cette protéine. Nous avons en outre pu montrer que chez *E. coli* également, le glycogène joue un rôle de stockage de sucre qui sert de source d'énergie pour faciliter la transition de la croissance bactérienne d'une source de carbone à une autre.

Le chapitre 7 s'intéresse à la dégradation de l'ensemble des ARNm d'E. coli. Je décrirai le développement d'un modèle simple reposant sur des approches de quasi-équilibre et permettant de prédire la cinétique de dégradation de chacun des ARNm cellulaires de la bactérie E. coli. Nous avons pu formuler de nouvelles hypothèses sur le rôle possible de la compétition entre ARNm pour leur fixation au dégradosome lors de l'adaptation de la croissance bactérienne à des changements environnementaux. Nous montrons également que ce mécanisme de compétition joue un rôle physiologique grâce à une approche de modélisation non linéaire à effets mixtes utilisant le modèle mécanistique de la dégradation des ARNm et des jeux de données de transcriptomique dynamique mesurant la cinétique de disparition des ARNm cellulaires.

Le chapitre 8 est dédié à des projets en cours et futurs sur l'analyse intégrative de la croissance microbienne et les stratégies d'allocation de resources des bactéries. Les travaux menés dans le cadre de ces projets vont définir mon activité scientifique dans les années à venir et celle de la future équipe-projet, en cours de création, dont je prendrai la direction.
Acknowledgements

First, I would like to thank Alexander Bockmayr, Matthieu Jules, and Vassily Hatzimanikatis, who kindly accepted to review this manuscript. Alexander guided my first steps in the modelling of biological systems a long time ago. I am delighted to have him on my habilitation committee. I would also like to thank the other members of my committee, Sandrine Charles, Muriel Cocaign-Bousquet, and Hans Geiselmann for their time and availability.

The present habilitation thesis is the result of the research that I have carried out within the Inria project-teams HELIX and IBIS. Many thanks go to my colleagues in these teams for warm and friendly atmosphere, and for the many scientific and non-scientific discussions. Special thanks go to François Rechenmann, Alain Viari, Hidde de Jong, Hans Geiselmann, Michel Page, Eugenio Cinquemani, Aline Marguet, and to the members of the BIOP team at LIPhy, in particular Corinne Pinel, and to Olivier Ali. I would especially like to thank all students and young researchers whom I have had the pleasure to (co-)supervise or otherwise collaborate with: Valentina Baldazzi, Sara Berthoumieux, Ismail Belgacem, Stefano Casagranda, Thibault Etienne, Jérôme Izard, Nils Giordano, Edith Grac, Manon Morin, Stéphane Pinhal, Pedro Tiago Monteiro, and Valentin Zulkower. I am also grateful to Inria for the support to carry out this pluridisciplinary research and for providing a great work environment.

The research presented in this manuscript is the result of a broad collaborative effort. I consider myself lucky to have worked with so many talented and nice people. In addition to the people I mentioned above, many others contributed to this research, directly or indirectly. I would like to thank, in particular, Jean-Luc Gouzé, Aline Métris, Jozsef Baranyi, Grégory Batt, Andreas Kremling, Tomas Gedeon, Laurent Trilling, Yohann Couté, and Myriam Ferro. For succinctness of the manuscript, I have had to omit some of the papers we have worked on together. Be sure, however, that I value our collaborations. Special thanks go to Jean-Luc Gouzé for our many scientific discussions along the years and for his kindness. Thanks to him, model reduction has become less mysterious. I would like also to thank Marie-France Sagot and all members from the Inria project-team ERABLE for our nascent collaboration on the analysis of metabolic graph models.

A significant part of my research is now carried out in tandem with Muriel Cocaign-Bousquet. I would like to warmly thank her for our scientific and non-scientific discussions. I value them so much! I am yearning to go back to Toulouse in the post-COVID times to pursue our brainstorming meetings on mRNA degradation and cellular metabolism, wrapped up at the end of the day by a "crêpe suzette". I am happy that you will be a member of the new project-team that we are creating. Thank you also to all TBI members, in particular to the BLADE team, for the friendly atmosphere and the many scientific exchanges.

And... I would never be where I am without my family and my friends... Warm and heartfelt thanks to you, especially to my parents and my sisters and brother! My deepest gratitude goes to Hidde, and to Quentin, Inès, and Siaka: you are illuminating my life.
<figure><img src="image_3.png" /></figure>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Context</td><td>1</td></tr><tr><td>1.2</td><td>My journey in systems biology</td><td>4</td></tr><tr><td>1.3</td><td>Overview</td><td>5</td></tr><tr><td>2</td><td>Model approximation and reduction</td><td>7</td></tr><tr><td>2.1</td><td>Deterministic modelling of biochemical network models</td><td>7</td></tr><tr><td>2.1.1</td><td>General form</td><td>7</td></tr><tr><td>2.1.2</td><td>Flux analysis</td><td>8</td></tr><tr><td>2.1.3</td><td>Kinetic models</td><td>9</td></tr><tr><td>2.2</td><td>Approximation of large kinetic models</td><td>12</td></tr><tr><td>2.2.1</td><td>Model reduction based on time-scale separation</td><td>12</td></tr><tr><td>2.2.2</td><td>Model approximation by means of piecewise-linear functions</td><td>16</td></tr><tr><td>2.2.3</td><td>Discussion and perspectives</td><td>18</td></tr><tr><td>3</td><td>Qualitative analysis of the dynamics of gene regulatory networks</td><td>21</td></tr><tr><td>3.1</td><td>Qualitative modelling and simulation of piecewise-linear models</td><td>21</td></tr><tr><td>3.2</td><td>Qualitative analysis of the carbon starvation response in E. coli</td><td>22</td></tr><tr><td>3.3</td><td>Discussion and perspectives</td><td>25</td></tr><tr><td>4</td><td>Analysis of dynamical gene expression and metabolomics data</td><td>29</td></tr><tr><td>4.1</td><td>Estimation of time-varying growth, uptake and secretion rates from dynamic metabolomics data</td><td>29</td></tr><tr><td>4.2</td><td>Estimation of promoter activities and protein concentration profiles from reporter gene data</td><td>32</td></tr><tr><td>4.3</td><td>Discussion and perspectives</td><td>36</td></tr><tr><td>5</td><td>Analysing and controlling cell physiology</td><td>39</td></tr><tr><td>5.1</td><td>Contribution of cell physiology to the global control of gene expression</td><td>39</td></tr><tr><td>5.2</td><td>A synthetic biology approach to control bacterial growth</td><td>42</td></tr><tr><td>5.3</td><td>Discussion and perspectives</td><td>45</td></tr><tr><td>6</td><td>Metabolic network models as platforms for integrating omics data</td><td>49</td></tr><tr><td>6.1</td><td>Post-transcriptional regulation of central carbon metabolism in E. coli</td><td>49</td></tr><tr><td>6.2</td><td>Post-transcriptional regulation of metabolic adaptation</td><td>52</td></tr><tr><td>6.3</td><td>Discussion and perspectives</td><td>56</td></tr><tr><td>7</td><td>Analysis of bacterial mRNA decay</td><td>57</td></tr><tr><td>7.1</td><td>Competitive effects in bacterial mRNA decay</td><td>57</td></tr><tr><td>7.2</td><td>Integrative analysis of mRNA degradation</td><td>60</td></tr><tr><td>7.3</td><td>Discussion and perspectives</td><td>62</td></tr></table>
# Introduction

## 1.1 Context

This manuscript is an overview of my research activities at Inria over the past fifteen years. The common denominator of the work presented is the use of computational systems biology approaches to unravel the complex molecular mechanisms involved in the adaptation of microorganisms to their environment. My work thus heavily relies on experimental data, seeking to maximally exploit and make sense of them. Before I describe my own contributions, I will place them in the context of the constantly evolving field with multiple facets that is systems biology. I will focus in particular on its historical context, because it has shaped the methodology and practice of modern systems biology in general and my own research trajectory in particular.

Much has been said and written about systems biology [Ideker et al., 2001, Kitano, 2001, 2002, Kirschner, 2005, Wolkenhauer, 2001, to name but a few]. Still trying to establish itself as a new field of science, systems biology has about as many definitions as there are systems biologists [Calvert and Fujimura, 2011]. We here summarize the literature by saying that systems biology refers to the body of approaches seeking to understand biological systems at the systems level. However, there is a general agreement on the agenda of systems biology, that it moves away from just cataloguing and characterizing cell components at the molecular level, to look at how they functionally interact and give rise to complex biological processes [Ideker et al., 2001, Kitano, 2001]. Reaching such understanding requires a collection of methods and ideas from various disciplines. It includes the biochemical and molecular biology approaches to analyse network components and their interactions. It shares with physiology the ultimate goal and methods to understand the functioning of living organisms and the underlying mechanisms, and with developmental biology, the analysis of how these mechanisms lead to a succession of physiological states. It makes use of various statistical approaches to analyse the vast amount of molecular data and of computational models of emergent phenomena for predicting the behaviour of biological systems. Eventually, systems biology includes aspects from evolutionary biology and ecology, acknowledging that living organisms are products of selection, which is more difficult to understand on a molecular level [Kirschner, 2005]. The combination of these various disciplines depends on the problem studied and allows to obtain an overall explanation of complex biological phenomena, which would not have been possible otherwise [Brigandt, 2010].

Systems biology is thus an integrative field that has been built on progresses in the disciplines mentioned above and past attempts. As such, it is both an old and a new field in biology [Auffray and Nottale, 2008, Kitano, 2002]. An old field, because system-level understanding in biology was already proposed between the 1930s and the 1970s. Far from being exhaustive, we can cite pioneering works such as the General System Theory introduced by the physiologist L. von Bertalanffy which, applied to biology, aimed at studying living organisms as a whole and considered them as open systems interacting with the environment [Rosen, 1958, von Bertalanffy,
1940, Von Bertalanffy, 1950, 1969]. We owe the notion of homeostasis to the physiologist W.B. Cannon [Cannon, 1929, 1939, Cooper, 2008], while the control engineer N. Wiener introduced the concept of negative feedback [Wiener, 1948], central to understand homeostasis. The name of systems biology itself was coined about fifty years ago [Mesarović, 1968] and first theoretical studies of morphogenesis [Turing, 1952], metabolism [Heinrich and Rapoport, 1974, Kacser, 1973, Mitchell, 1961, Rottenberg et al., 1967, Savageau, 1976], gene expression [Goodwin et al., 1963, Kauffman, 1971], and the nonlinearities of complex biological systems [e.g., Prigogine and Nicolis, 1971] were published over the same period.

All these studies laid the ground for the systems view in biology, but systems biology itself took off only at the beginning of the twenty-first century as a new biological field [Auffray and Nottale, 2008, Kitano, 2002]. Various explanations have been offered. First, early works like those of L. von Bertalanffy were considered too abstract and biological systems too complex to be modelled. As an illustrative example of the latter point, I will cite the French physiologist Claude Bernard. In the nineteenth century, he already advocated the use of mathematics in biology [Bernard, 1865, p238]:

> "Cette application des mathématiques aux phénomènes naturels est le but de toute science, parce que l'expression de la loi des phénomènes doit toujours être mathématique."

However, the idea was already there that the complexity of biological phenomena and incomplete knowledge prevent their mathematical study [Bernard, 1865, p238-239]:

> "Or je pense que les tentatives de ce genre sont prématurées dans la plupart des phénomènes de la vie, précisément parce que ces phénomènes sont tellement complexes, qu'à côté de quelques-unes de leurs conditions que nous connaissons, nous devons non seulement supposer, mais être certain, qu'il en existe une foule d'autres qui nous sont encore absolument inconnues. [...] Ce n'est point que je condamne l'application mathématique dans les phénomènes biologiques, car c'est par elle seule que, dans la suite, la science se constituera; seulement j'ai la conviction que l'équation générale est impossible pour le moment, l'étude qualitative des phénomènes devant nécessairement précéder leur étude quantitative."

Other explanations offered for the belated emergence of systems biology are essentially a matter of agenda. Somehow, they have followed the idea of Claude Bernard that an accumulation of detailed information at the physiological, cellular, and molecular levels is needed beforehand. An encyclopedic knowledge of the genes and proteins involved in biological phenomena has thus been developed over the years. On a small scale and through *ad hoc* studies at the beginning, but things have accelerated with the availability of genome sequences and the development of high-throughput approaches [Ideker et al., 2001, Kitano, 2002, Voit, 2017]. They allowed the inventory and quantification of cell components, in a wave of technological developments improving the probing, sensing, and imaging of biological systems. Parallel developments in mathematics, physics, and computer science have made theoretical and simulation tools more powerful and more accessible to a wide audience. Hence, following the development of bioinformatics for the processing and statistical analysis of the vast amount of high-throughput data, systems biology has allowed to make sense of the data, by analysing in
terms of biological networks the relationships between genes, proteins, and their function.

I took the time to survey the history of systems biology to better illustrate why systems biology as we know it nowadays is built on the two pillars of experimental data analysis and mathematical modelling, the so-called top-down and bottom-up approaches. The two approaches are the fruit of the evolution of the field and at the centre of an epistemological debate [Ideker et al., 2001, Kell and Oliver, 2004, Westerhoff and Palsson, 2004]. Top-down systems biology is more related to the modern branch of the domain, with the expansion of molecular biology to genome-wide analyses. A data-driven approach in essence, it is concerned with the analysis of large-scale data sets, whose access has been greatly facilitated by the recent technological, algorithmic and computational advances [Bruggeman and Westerhoff, 2007, Ideker et al., 2001]. The challenge is to gain a system-level understanding by integrating these omics data. The induction task is not trivial due to the fact that data can be multi-variate and multi-layered (e.g. transcriptomics, metabolomics, and proteomics data). The most commonly used models to that aim are often phenomenological, in the sense that they neither describe regulatory mechanisms nor include prior knowledge of the molecular components. They account for correlations between measured concentrations, which allows identifying groups of co-regulated genes and formulating hypotheses on biological mechanisms at work, although correlation does not necessarily mean causality [Bersanelli et al., 2016, Noor et al., 2019]. Top-down systems biology is thus more observational, close to the "naturalist" attitude, and – for its detractors – a "fishing expedition" [Calvert and Fujimura, 2011, Kell and Oliver, 2004].

At the other side of the spectrum, the bottom-up approach is hypothesis driven and more related to the old débuts of the field in the mid-twenties century. It deduces the functional properties of a subsystem, through the development of mathematical models, which describe mechanistically how components interact and predict how the system behaviour emerges from these interactions [Bruggeman and Westerhoff, 2007]. This approach embraces the aspirations of physics and engineering, with the idea to uncover laws and make the behaviour of biological systems predictable and controllable [Calvert and Fujimura, 2011, Kell and Oliver, 2004, Westerhoff and Palsson, 2004]. This has paved the way for the design of biological systems that function close to specification, at the centre of the synthetic biology agenda [Arkin, 2013, Smolke and Silver, 2011].

While top-down and bottom-up approaches have been often opposed, they are just two sides of the same coin, with the common goal of relating the phenotype to the genotype [Kell and Oliver, 2004, Westerhoff and Palsson, 2004]. The literature in the field reflects a progressive reconciliation of these approaches. Constraint-based modelling of metabolic networks is one such example. A favourite approach to model metabolic networks, it provides a comprehensive representation of the metabolism of an organism, by representing the metabolic network through a series of physico-chemical constraints, including reaction stoichiometry and assumption of steady-state metabolite concentrations [Price et al., 2004, Orth et al., 2010, Volkova et al., 2020]. Even though it requires prior knowledge of the network connectivity, constraint-based modelling pertains more to a top-down approach in my opinion, because it makes extensive use of genomic, transcriptomic, proteomic, and metabolomic data and aids their interpretation [Lewis et al., 2010, Machado and Herrgård, 2014, Ramon et al., 2018, Shlomi et al., 2008, Volkova et al., 2020]. However, constraint-based modelling in its classical implementation is unable to give an
insight into cellular substrate concentrations. That is what kinetic modelling does, but this approach suffers because parametrizing mechanistic models is both costly and time-consuming. Bridging the gap between constraint-based models and kinetic models, and thus between top-down and bottom-up approaches, has been a natural expansion of these approaches. Subsequent improvements to constraint-based methods have allowed us to quantitatively understand metabolism and its regulation, notably by accounting for metabolite concentrations [Henry et al., 2007, Hoppe et al., 2007, Kümmel et al., 2006a] and coupling constraint-based models to mechanistic models [e.g., Cotten and Reed, 2013, Covert et al., 2008, Hanly and Henson, 2011, Lee et al., 2008, Mahadevan et al., 2002, Smallbone et al., 2007, Yizhak et al., 2010]. These models predict dynamic metabolic behaviour without explicitly solving ordinary differential equations. Most recent developments allow the parametrization of constraint-based models, turning them into nearly genome-scale kinetic models [Gopalakrishnan et al., 2020].

The research activities presented in this manuscript reflect these different trends of systems biology, with work relying on the use of mathematical models describing the dynamics of biochemical networks and other works related to the model-based interpretation of high-throughput data to explore and characterize molecular mechanisms involved in bacterial adaptation to environmental changes.

## 1.2 My journey in systems biology

Trained as a biochemist and molecular biologist at the University of Nancy, my scientific interest has always been in interdisciplinary studies of the functioning of biological systems. Over the course of the years, my research has gradually evolved from experimental biology to computational systems biology, although at times I continue to run wet-lab experiments to answer specific research questions or to acquire new data. Throughout this journey, a major concern has been to obtain good data and/or to make it informative for a tighter integration of models and experiments. The nature of the data, its quality and quantity, were determinant for the methodological choices in all studies. In this sense my work has followed the evolution of systems biology and that of biology becoming more and more quantitative.

My first steps in research were in enzymology, during an internship for my Bachelor's degree. I studied the allosteric regulation of a glycolytic enzyme (GAPDH) of *Bacillus stearothermophilus* using a combination of molecular biology and biochemistry tools. As a master and PhD student, I focused on the regulation of HIV-1 RNA splicing by human proteins. HIV-1 expresses its gene products from a single transcript, which undergoes alternative splicing using a combination of four donor and nine acceptor sites. By means of molecular biology, biochemistry, and bioinformatics tools, I studied the regulatory mechanisms driving the choice of acceptor sites by the host cell splicing machinery [Jacquenet et al., 2001, Ropers et al., 2004, Hallay et al., 2006, Khoury et al., 2009]. In collaboration with Alexander Bockmayer and Damien Eveillard, then at Inria in Nancy, I could analyse by means of mathematical modelling and constraint-programming tools how the competition between activating and inhibitory proteins for their binding in the vicinity of splicing sites determines their usage by the splicing machinery [Eveillard et al., 2003, 2004]. The difficulty consisted in the complexity of RNA splicing and the lack of quantitative information to model the process. This first exposure to mathematical modelling shaped the years to come. I pursued the experience by joining the Helix Bioinformatics group at Inria
Grenoble - Rhône-Alpes as a post-doctoral researcher. In collaboration with Hidde de Jong and Hans Geiselmann (now with the interdisciplinary laboratory of physics, LIPhy), I used qualitative modelling approaches to study the functioning of a genetic regulatory network controlling the carbon starvation response of the bacterium *Escherichia coli*. The mathematical formalism allowed to compensate for the lack of quantitative information, although stress responses were relatively well characterized in *E. coli*. The following years witnessed a rapid growth of dynamical and quantitative data at all cell levels, which opened new avenues for integrative analyses of bacterial adaptation to environmental or genetic cues, looking not only at the gene expression level but also at metabolism and cell physiology. This motivated me to adopt more quantitative modelling approaches to analyse the relation between genotype and phenotype. Now as a research scientist with the systems biology group Ibis, my research activities cover aspects of data analysis, mathematical modelling of biochemical networks, and cell growth, with a tight integration of heterogeneous and multi-layered experimental data.

Analysing complex biochemical systems by means of models and experiments requires a broad expertise. On the biological side, the most important skills are in biochemistry, molecular biology, and microbiology. A good knowledge of the biological systems and of the experimental methods is required for mathematical modelling. On the theoretical side, the major domains of expertise are biostatistics, bioinformatics, and dynamical systems. These are critical for data analysis, model development, analysis, and identification. One of my main objective was to gain know-how in these various domains to achieve my research agenda. I enlarged my spectrum of expertise through various collaborations. Currently, part of my research is nurtured by challenging problems on the regulation of cellular metabolism and mRNA metabolism in particular, in the context of an active collaboration with Muriel Cocaign-Bousquet and her group at the Toulouse Biotechnology Institute. The biochemical models of gene expression and cell growth that I develop quickly grow in complexity and require more advanced tools for their reduction and dynamical analysis. I work on these aspects in collaboration with Jean-Luc Gouzé at Inria Sophia Antipolis - Méditerranée. All these models rely on the integration of data. Their analysis and use in parameter estimation pose challenging statistical problems that I tackle in collaboration with Eugenio Cinquemani within Ibis. The experimental aspects of my research, the synthetic biology application, and the modelling of gene expression and cell growth thrive on many interactions with Hidde de Jong and Hans Geiselmann over the course of the years.

1.3 Overview

In the remainder of this manuscript, I will present a selection of my research activities since
I joined Inria. I will focus on bacterial growth adaptation, omitting the works dedicated to
HIV-1 RNA splicing published after my PhD thesis [Eveillard et al., 2004, Ropers et al., 2004,
Hallay et al., 2006, Khoury et al., 2009]. In Chapter 2, I will discuss mathematical models
to investigate the dynamical functioning of gene regulatory networks, first where quantitative
data are poor if not absent. Models in this case are qualitative, resulting from simplifications
of more complex algebro-differential models that can be used to quantitatively describe better-
characterized systems. The latter models typically combine different time scales, which poses
interesting questions that we addressed by model reduction and simplification to obtain models
that are a good compromise between simplicity and biological realism. Chapter 3 is dedicated
to the application of qualitative models of gene expression to analyse the carbon starvation
response of *Escherichia coli*. The development of biology into a quantitative science requires the
development of approaches for the analysis of metabolomics data as well as reporter gene data.
This will be the topic of Chapter 4. In Chapter 5, I will discuss work related to the analysis and
control of cell physiology, and how the study of global mechanisms of control of gene expression
led to the development of a synthetic strain of *E. coli* allowing to shift the allocation of resources
from growth to production of a metabolite of interest. Chapter 6 is dedicated to the analysis
of metabolic networks. Metabolomics and other low to high-throughput data are more and
more accessible nowadays, which makes metabolic network models suitable platforms for the
integration of experimental data to analyse the regulation of metabolism in various genetic and
environmental backgrounds. Chapter 7 presents an analysis of bacterial mRNA decay at the
genome-wide level, taking advantage of the availability of dynamical transcriptomic data. The
outlook in Chapter 8 is the opportunity to conclude and discuss the future. The annex is
concerned with administrative and vitae information, as well as a list of the papers highlighted
in the manuscript with their abstract.
# Model approximation and reduction

Mathematical models must be realistic representations of the biological system, from the general topology of the underlying biochemical network to the regulation exerted on the biological processes. However, the modelling task is far from trivial as knowledge of the network topology is not exhaustive and quantitative information is often missing. Different modelling formalisms exist to overcome these limitations and adjust the level of precision to the available data. Here we will restrict the description to deterministic approaches. They are well adapted to study network dynamics in cell populations, which is central to my research activities.

The first section introduces basic concepts of biochemical network modelling. I will notably focus on flux analysis approaches, which will be central in Chapter 6. The description of kinetic modelling approaches will develop approximations used in Michaelis-Menten kinetics, their domain of validity, and issues related to parameter estimation. The analysis of mRNA decay described in Chapter 3 relies on these concepts. In the second section, I will describe some of my work on the approximation and reduction of biochemical models, based on a paper in *IEEE/ACM Transactions on Computational Biology and Bioinformatics* [Ropers et al., 2011]. The approaches described in this article were repeatedly used in the modelling works described elsewhere in the manuscript.

## 2.1 Deterministic modelling of biochemical network models

### 2.1.1 General form

Ordinary differential equations are the classical formalism for modelling the dynamics of natural and man-made systems. In the case of biochemical reaction networks, they relate the rate of change of variables to their values by means of mass-balance equations. This gives in matrix notation [Heinrich and Schuster, 1996]:

$$ \frac{dx}{dt} = Nv - \mu x, \quad x(0) = x_0. \qquad (2.1) $$

In this system of coupled equations, $x \in \mathbb{R}_+^n$ and $v : \mathbb{R}_+^n \to \mathbb{R}^q$ and denote the vectors of reaction rates and concentrations at time $t$, respectively, $N \in \mathbb{Z}_+^{n \times q}$ is the stoichiometry matrix, and $\mu \in \mathbb{R}_+$ the growth rate. This type of model is based on the premise that the system is well-stirred and that concentrations can be regarded as continuous quantities. We consider here constant pH and temperature conditions because these two parameters are generally controlled in the experiments. The model is used to describe various biological processes such as transcription, translation, enzymatic reactions, complex formation, and transport reactions between compartments.
2.1.2 Flux analysis

When the biological system under study is at steady state, the balance equation in (2.1) becomes a system of algebraic equations [Heinrich and Schuster, 1996]:

$$Nv - \mu x = 0. \qquad (2.2)$$

By convention, reactions rates at steady-state are termed fluxes. This mathematical description is typically used for metabolic networks, in which case the flux vector $v$ and the stoichiometry matrix $N$ are restricted to the internal metabolites.

An additional and frequently made assumption is that growth dilution of intracellular metabolites is negligible with respect to the turn-over of metabolite pools by enzymatic reactions:

$$Nv = 0. \qquad (2.3)$$

The system (2.3) is commonly used to analyse stationary fluxes. There is no longer an explicit dependence of fluxes on concentrations and the fluxes are the new variables of the system.

While we may have measurements for some of them, the vast majority of fluxes in vector $v$ is unknown and makes the system under-determined. It is solved by means of constraint-based modelling approaches. Basically, equation (2.3) constrains fluxes to the null space of the stoichiometry matrix $N$. There is generally no unique solution to this system but a distribution of feasible stationary solutions. It can be narrowed down by the addition of constraints in the form of inequalities applied to the fluxes of each reaction $i$ [for review, Bordbar et al., 2014]:

$$lb_i \leq v_i \leq ub_i. \qquad (2.4)$$

For instance, the range of incoming and outgoing fluxes can be fixed on the basis of measured uptake and secretion rates [Mo et al., 2009]. Genetic knock-outs are carried out by setting the bounds of the associated reactions to zero [Edwards et al., 2001]. Thermodynamic constraints allow to eliminate infeasible cycles or to set reaction directions based on measured intracellular metabolite concentrations and the Gibbs free energy of formation [Beard et al., 2004, Henry et al., 2007, Hoppe et al., 2007, Kümmel et al., 2006b, Müller and Bockmayr, 2013, Qian and Beard, 2005]. These simple constraints have evolved to more advanced approaches with the growing availability of omics data. Transcriptomics and proteomics data can be used to fix reaction bounds to zero for reactions corresponding to absent mRNAs or proteins or by linearly adjusting the bounds to the mRNA or protein abundances [Åkesson et al., 2004, Chandrasekaran and Price, 2010, Colijn et al., 2009, Covert et al., 2001, Tian and Reed, 2018, van Berlo et al., 2009, Yizhak et al., 2010]. The use of omics data allows the development of context-specific models [e.g. Agren et al., 2012, Bordbar et al., 2014, Heirendt et al., 2019, J senior et al., 2020, Jensen and Papin, 2011, Shlomi et al., 2008, Thiele et al., 2020, Wang et al., 2012]. This improves the quality of the predictions and turn metabolic network models into scaffolds for the analysis of high-throughput data. Available information on kinetic parameters can also be used to tighten bounds on reaction fluxes [Cotten and Reed, 2013, Fleming et al., 2010].

The application of constraints to reaction fluxes shrinks the solution space to a biologically relevant region. Possible steady-state flux distributions within the region form a convex
polyhedral cone in a high-dimensional space that can be analysed by different but related approaches. Methods for the analysis of pathways, based on elementary mode analysis [Schuster et al., 1999] or extreme-pathway analysis [Schilling et al., 2000], are used to define the limitations and production capabilities of metabolic systems. Flux balance analysis (FBA) uses linear programming to select an optimal flux distribution within the region, for instance by maximizing an objective function such as the biomass or ATP production [Bonarius et al., 1997, Edwards and Palsson, 1999, Sauer et al., 1998, Varma and Palsson, 1994a,b]. The main difference with metabolic flux analysis is the use in FBA of an objective function [Wiechert, 2001, Zupke and Stephanopoulos, 1994]. Flux variability analysis (FVA) has been developed to identify the possible alternative optima of the FBA solution [Mahadevan and Schilling, 2003]. It determines the maximum and minimum fluxes through each reaction when the flux of the objective function is constrained to its maximum value. The choice of an optimization function in FBA and FVA is not trivial, although optimizing cell growth and energy use has been shown to correctly predict metabolic fluxes in microorganisms [Carlson and Srienc, 2004a,b, Edwards and Palsson, 2000, Edwards et al., 2001, Feist and Palsson, 2010, Ibarra et al., 2002]. Studies have questioned the universality of the assumption of such a metabolic objective [Harcombe et al., 2013, Molenaar et al., 2009, Schuetz et al., 2007]. An unbiased alternative approach is random sampling. It explores the solution space with Monte Carlo approaches without optimizing an objective and returns probability distributions of the fluxes [Haraldsdóttir et al., 2017, Kaufman and Smith, 1998, Keaty and Jensen, 2020, Megchelenbrink et al., 2014, Wiback et al., 2004]. The sampling approach is more and more applied for the analysis of metabolic networks [Bordbar et al., 2014].

Metabolic network modelling is an active field of research with an ever growing range of applications and methodologies. In this section, I only sketched the main constraint-based modelling approaches for the sake of brevity. These approaches make use of genome-scale models of cell metabolism, whose reconstruction is far from trivial and time-consuming [Gu et al., 2019, Gudmundsson et al., 2017]. Fortunately, numerous reconstructions are freely available in databases such as BIGG [King et al., 2016, Norsigian et al., 2020]. The models are encoded and annotated in the standardized SBML format that facilitates their exchange and simulation with the above-mentioned methods on various platforms such as the COBRA toolbox [Becker et al., 2007, Heirendt et al., 2019]. We will use the latter for our analyses of cell metabolism in Chapter 6.

### 2.1.3 Kinetic models

Beyond metabolic network modelling, stoichiometric models can describe a whole range of different biochemical networks using the form (2.1). Taking into account the interactions between genes, proteins, and RNAs requires one to consider the kinetics of the reactions and their dependence on molecular concentrations and kinetic parameters:

$$ \frac{dx(t)}{dt} = Nv(x, p) - \mu x, \quad x(0) = x_0, \qquad (2.5) $$

with $p \in \mathbb{R}_+^k$ a vector of parameters. The vector of rate laws $v(x, p)$ is a generally nonlinear function of concentrations and parameters.
A variety of mathematical functions exist for kinetic rate laws. The particular choice of a mathematical form depends on the degree of knowledge that we have on the system and the level of precision that we wish. In the absence of precise knowledge, a standard rate law can be chosen. The law of mass action is one example [for a review, see Voit et al., 2015]. It describes the rate $v_j$ of an elementary chemical reaction $j$ (with a single mechanistic step) as being proportional to the product of the concentrations $x_i$ with a given stoichiometry $a_i$:

$$v_j = k_j \prod_i x_i^{a_i}, \qquad (2.6)$$

where $k_j$ is called the rate constant of the reaction. It states that increasing the number of molecules of a given species will increase the probability that they collide. The main advantage of models based on the law of mass action is that they can be determined directly from the elementary reactions and their stoichiometry. Enzymatic reactions are one example. For instance, a free enzyme E catalyses the transformation of a substrate S into a product P by first binding to the free substrate to form an enzyme-substrate complex C. In the example below, the latter can be converted irreversibly into the original enzyme and the product or, if the transformation fails, into the original enzyme and substrate:

$$E + S \underset{k_{-}}{\overset{k_{+}}{\rightleftharpoons}} C \overset{k_{cat}}{\rightarrow} E + P, \qquad (2.7)$$

where $k_+$ and $k_-$ are the reaction rate constants and $k_{cat}$ the catalytic constant of the reaction. Applying the mass-action law to the reaction in (2.7), we write the following system of differential equations

$$\begin{aligned} \frac{dx_E}{dt} &= -k_+ x_E x_S + (k_- + k_{cat}) C, \\ \frac{dx_S}{dt} &= -k_+ x_E x_S + k_- C, \\ \frac{dx_C}{dt} &= k_+ x_E x_S - (k_- + k_{cat}) C, \\ \frac{dx_P}{dt} &= k_{cat} x_C, \end{aligned} \qquad (2.8)$$

and initial conditions: $x_E(0) = x_S(0) = x_C(0) = 0, x_P(0) = 0$. Here it is assumed that no product and complex are present at the start. The total concentration of enzyme and substrate is conserved along the reaction: $x_E^0 = x_E + x_C$ and $x_S^0 = x_S + x_C + x_P$.

A popular model derived from these equations by means of the standard quasi-steady-state approximation (sQSSA) is the Henri-Michaelis-Menten equation [Briggs and Haldane, 1925, Cornish-Bowden, 2015, Michaelis and Menten, 1913]. The sQSSA assumes that the concentration of complex C rapidly equilibrates to its quasi-steady-state value. By solving $dx_C/dt = 0$ in (2.8), and using the mass-conservation relation for the enzyme concentration, we obtain the popular form of the equation describing the rate of accumulation of the product with saturation effects:

$$\frac{dx_P}{dt} = \frac{k_{cat} x_E^0 x_S}{K_m + x_S}, \qquad (2.9)$$

where $K_m = \frac{k_- + k_{cat}}{k_+}$ is the Michaelis-Menten constant.

The Michaelis-Menten equation has become a canonical approach to understand enzyme kinetics. It is widely used to estimate $k_{cat}$ and $K_m$ values from product progress curves [Duggleby and Clarke, 1991, Johnson, 2013, Stroberg and Schnell, 2016, Tummler et al., 2014]. The
in-vitro experiments are carried out in conditions of validity of the equation, when the total substrate concentration is much larger than the total enzyme concentration, so that the amount of substrate bound to the enzyme is negligible. This condition, which has been generalized by Segel [1988] and Segel and Slemrod [1989], implies that the quasi-steady-state assumption for the enzyme-substrate complex is valid at a lower enzyme concentration with respect to the sum of the total substrate concentration and $K_m$ value [for review, Schnell and Maini, 2003]:

$$ \frac{x_E^0}{K_m + x_S^0} \ll 1. \qquad (2.10) $$

Validity of the Michaelis-Menten equation is a necessary but not sufficient condition for an accurate estimation of kinetic parameter values [Chen et al., 2010, Choi et al., 2017, Stroberg and Schnell, 2016]. The highly correlated structure and non identifiability of the parameters requires a proper design of dynamic experiments, such as choosing appropriate times for data collection and an initial substrate concentration equal to (as a rule of thumb) two or three times the $K_m$ value [Choi et al., 2017, Duggleby and Clarke, 1991, Stroberg and Schnell, 2016]. The whole difficulty is that a proper assay thus requires a priori knowledge of the $K_m$ value... to determine the $K_m$ value more precisely.

An alternative approximation has been proposed for system (2.8). It overcomes the parameter estimation issues [Choi et al., 2017] and the observation that, in vivo, enzyme concentrations are often higher than in enzymatic assays or at least of the same magnitude as their substrate [e.g. Albe et al., 1990]. It is based on the total quasi-steady-state approximation (tQSSA) by replacing the free substrate concentration with its total substrate concentration. The derivation gives a somewhat complex expression [Borghans et al., 1996, Tzafriri, 2003]:

$$ \frac{dx_P}{dt} = k_{cat} \frac{x_E^0 + K_m + x_S^0 - x_P}{2} - \frac{\sqrt{(x_E^0 + K_m + x_S^0 - x_P)^2 - 4 x_E^0 (x_S^0 - x_P)}}{2}. \qquad (2.11) $$

From expression (2.11), Tzafriri [2003] developed a simpler approximation called first-order tQSSA. It resembles the sQSSA form of Michaelis-Menten kinetics:

$$ \frac{dx_P}{dt} = \frac{k_{cat} x_E^0 x_S}{K_m + x_E^0 + x_S}, \qquad (2.12) $$

and has the following criteria of validity:

$$ x_E^0 + K_m \gg x_S^0 \quad \text{and} \quad K_s \ll K_m, \quad \text{or:} \qquad (2.13) $$

$$ x_E^0 \gg x_S^0 \quad \text{and} \quad x_E^0 \gg K_m \approx K_s, \qquad (2.14) $$

where $K_s$ is the Van Slyke-Cullen constant defined by $K_s = k_{cat}/k^+$. Compared to sQSSA, the tQSSA approximation extends the parameter domain for which it is reasonable to assume that $dx_C/dt \approx 0$ [Schnell and Maini, 2000, 2003, Tzafriri, 2003]. With this approximation, more accurate and precise estimations can be obtained with a proper experimental design and without prior information on parameter values [Choi et al., 2017]. The tQSSA approximation has been used for modelling different types of biochemical systems [e.g. Ciliberto et al., 2007,
Pedersen et al., 2008b,a, Tzafiri et al., 2002]. We also used it in a study of mRNA degradation kinetics described in Chapter 7.

Various other quasi-steady-state approximations have been proposed, as well as different ex-
tensions of those described here [Schnell and Maini, 2003]. In addition, many other expressions
allow the modelling of the kinetics of enzymatic reactions, with inhibition, allostery,... or of gene
expression with Hill functions for instance. It is beyond the scope of this section to describe
these in detail, as there are excellent textbooks on this topic [e.g. Klipp et al., 2016, Segel,
1993, Voit, 2017]. The important message is that the mass-action law, the Michaelis-Menten
equation and similar rate laws offer a framework to model a large diversity of biochemical
processes, from enzymatic reactions to gene expression. This allows to obtain ODE models for
the solution of which numerous numerical methods and software tools are available.

ODE models of biochemical systems are nonlinear, often stiff as they include processes evolv-
ing on different time scales, and their size grows quickly as more biochemical processes are con-
sidered. As a consequence, they are not amenable to formal mathematical analysis and reliable
parameter estimation is a difficult problem due to a generalized lack of data. As we have seen in
this section, the problem can be partially solved with a proper experimental design and by using
a combination of mass conservation relations and valid quasi-steady-state approximations to
reduce the model dimension. In the following section, I will present some personal work related
to model approximations of larger biochemical systems and their validity, with two different sets
of approaches: model reduction based on time-scale separation and the use of quasi-steady-state
approximations, and model approximation by piecewise-linear functions, which is useful when
no quantitative data is available.

## 2.2 Approximation of large kinetic models

Various approximations have been proposed in the literature to reduce the size and complexity
of biochemical network models, tailored to typical response functions and time-scale hierarchies
found in genetic or metabolic regulation [de Jong and Ropers, 2006b, Heijnen, 2005, Heinrich
and Schuster, 1996, Okino and Mavrovouniotis, 1998, Papin et al., 2004, Pecou, 2005, Radulescu
et al., 2012, Roussel and Fraser, 2001, Savageau, 2001]. The approximations result in models
that are easier to handle mathematically and computationally, while maintaining important
properties of the original system. In particular, they reduce the dimension and the number of
parameters, and simplify the mathematical form of the equations.

2.2.1 Model reduction based on time-scale separation

This section is concerned with model-order reduction. The principle is based on the identification
of relationships between model variables such that fewer species need to be measured exper-
imentally or simulated. The derivation of the Michaelis-Menten equation is one prototypical
example: starting from the mass-action model with four variables, we arrive after reduction at a
single ODE describing the time evolution of the product concentration. Product formation can
be monitored experimentally, but neither the concentration of the enzyme-substrate complex nor
the concentrations of free substrate and free enzyme can be measured. The problem is that many
biochemical systems of interest combine a variety of processes, not only enzymatic reactions, but also gene expression, signalling... The reduction of such large and stiff systems is all but intuitive.

Three general strategies haven been pursued for model-order reduction [Okino and Mavrovouniotis, 1998]: (i) lumping, which transforms the original variables into a vector of lower dimension; (ii) sensitivity analysis neglecting network reactions and species with small impact; and (iii) time-scale analysis identifying the different time scales on which the network species react and considering the fast time-scale reactions and species at steady state. The latter approach is widely used in enzyme kinetics, as we have seen before, and more generally in biochemical network modelling with the usage of quasi-steady-state approximations. The approach has been mathematically formalized by singular perturbation theory [O'Malley, 1991]. We will focus on this approach in what follows. The objective is to obtain models of smaller dimension, with less variables and parameters, while the remaining variables and parameters are easier to measure experimentally. In addition, the parameters must be physiologically meaningful.

Given the ODE model as in (2.5), the first reduction step is to identify the different time scales of the system. This requires *a priori* knowledge such as the values of the parameters or other information about the biochemical processes. For instance, bacterial pools of metabolites have a turn-over of the order on second, minutes for mRNA pools, and hours for proteins [Shamir et al., 2016]. This is sufficient information in numerous cases. Consider for instance the small biochemical network in Figure 2.1(a). It is part of larger networks that we modelled and studied in [Baldazzi et al., 2010, 2012, Ropers et al., 2006, 2011]. The depletion of a carbon source like glucose is signalled to the cell by the phosphotransferase system. It leads to the activation of the adenylate cyclase Cya and the production of cAMP by the enzyme. This signalling molecule makes a complex with the transcription factor CRP, which binds to promoter regions and activates or inhibits transcription. This small module, which we will later refer to as an activation network, allows the gene expression program to be remodelled so that cells can cope with carbon starvation. In this small network, we identify fast processes, cAMP production and complex formation between CRP and cAMP. The slow processes correspond to the synthesis of CRP and Cya. The corresponding ODE model is shown in Figure 2.1(b). It includes five variables, describing the free concentrations of CRP, Cya, and cAMP, and the concentration of bound Cya and CRP. Rate laws include the Michaelis-Menten equation for cAMP synthesis, the mass-action law for the rate of association and dissociation of CRP·cAMP, while Hill functions are used to represent the regulation of CRP and Cya synthesis by the complex. Its general form is similar to (2.5).

If we are generally able to identify which processes are fast and which ones are slow, difficulties may arise for the separation of variables into fast and slow. We encounter this situation with the ODE model of the activation network shown in Figure 2.1(b). For instance the state equation for the free CRP concentration depends on fast processes, the association and dissociation of the complex CRP·cAMP, and slow processes such as the synthesis of the protein and its degradation. A linear transformation of the variables is needed to uncover the two time scales [Heinrich and Schuster, 1996]. We introduce vectors of slow and fast variables, $x^s \in \mathbb{R}_+^m$ and $x^f \in \mathbb{R}_+^{n-m}$,
<figure><img src="image_4.png" /></figure>

$$ \dot{x}_{y\sim} = \kappa_y^1 + \kappa_y^2 h^+(x_{c\sim m}, \theta_{c\sim m}^2, m_{c\sim m}) - \gamma_y x_{y\sim} + (k_{-1} + k_2 h^+(u_s, \theta_s, m_s)) x_{y\sim p} - k_1 x_{y\sim p} $$

$$ \dot{x}_{c\sim} = \kappa_c^1 + \kappa_c^2 h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m}) - \gamma_c x_{c\sim} + k_{-4} x_{c\sim m} - k_4 x_{c\sim} x_{m\sim} $$

$$ \dot{x}_{y\sim p} = k_1 x_{y\sim} x_p - (k_{-1} + k_2 h^+(u_s, \theta_s, m_s) + \gamma_y) x_{y\sim p} $$

$$ \dot{x}_{c\sim m} = k_4 x_{c\sim} x_{m\sim} - (k_{-4} + \gamma_c) x_{c\sim m} $$

$$ \dot{x}_{m\sim} = k_2 h^+(u_s, \theta_s, m_s) x_{y\sim p} + k_{-4} x_{c\sim m} - k_3 x_{m\sim} - k_4 x_{c\sim} x_{m\sim} $$

(b)

Slow system

$$ \dot{x}_y = \kappa_y^1 + \kappa_y^2 h^+(x_{c\sim m}, \theta_{c\sim m}^2, m_{c\sim m}) - \gamma_y x_y $$

$$ \dot{x}_c = \kappa_c^1 + \kappa_c^2 h^-(x_f, \theta_f^2, m_f^2) h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m}^1) + \kappa_c^3 h^-(x_f, \theta_f^1, m_f^1) - \gamma_c x_c $$

Fast system

$$ \dot{x}_{y\sim p} = k_1 (x_y - x_{y\sim p}) x_p - (k_{-1} + k_2 h^+(u_s, \theta_s, m_s) + \gamma_y) x_{y\sim p} $$

$$ \dot{x}_{c\sim m} = k_4 (x_c - x_{c\sim m}) (x_m - x_{c\sim m}) - (k_{-4} + \gamma_c) x_{c\sim m} $$

$$ \dot{x}_m = k_2 h^+(u_s, \theta_s, m_s) x_{y\sim p} + (k_3 - \gamma_c) x_{c\sim m} - k_3 x_m $$

(a)

<figure><img src="image_5.png" /></figure>

(c)

$$ \dot{x}_y = \kappa_y^1 + \kappa_y^2 - \gamma_y x_y $$

$$ \dot{x}_c = \kappa_c^1 + \kappa_c^2 h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m}) - \gamma_c x_c $$

$$ x_{c\sim m} = \frac{h^+(u_s, \theta_s, m_s) x_c x_y}{K_4 K_3 + h^+(u_s, \theta_s, m_s) x_y} $$

$$ K_3 = \frac{k_3}{k_2} $$

$$ K_4 = \frac{k_{-4}}{k_4} $$

(d)

(e)

**Figure 2.1** – (a) Activation network (from Ropers et al., 2011). (b) Detailed ODE model for the activation network. $x_{y\sim}$, $x_{c\sim}$, $x_{y\sim p}$, $x_{c\sim}$, and $x_{m\sim}$ denote the concentrations of free Cya, free CRP, Cya-ATP, CRP-cAMP, and free cAMP, respectively, while $u_s$ denotes the external glucose concentration. The total concentrations of Cya, CRP, and cAMP are referred to as $x_y$, $x_c$, and $x_m$, respectively. $h^+$ denotes a positive Hill function: $h^+(x, \theta, m) = x^m/(x^m + \theta^m)$. (c) Time-scale separation. (d) Model solutions. The blue curve represents a solution for the concentration variable $x_{c\sim m}$ in the original model. The red curve is the corresponding solution for $x_{c\sim m}$ in the QSS model. After an initial transient the solution of the original model rapidly relaxes to the QSS solution. (e) QSS model for the activation network. The model approximates the original model by coupling the fast variable $x_{c\sim m}$ to the slow variables $x_y$ and $x_c$. The slow variables are defined as: $x_y = x_{y\sim} + x_{c\sim m}$; $x_c = x_{c\sim} + x_{c\sim m}$; $x_m = x_{m\sim} + x_{c\sim m}$.
respectively ($m < n$). These are defined as linear combinations of the original variables $x$:

$$ \begin{bmatrix} x^s \\ x^f \end{bmatrix} = T x, \qquad (2.15) $$

with $T \in \mathbb{Z}^n \times \mathbb{Z}^n$. This gives for the stoichiometry matrix:

$$ \begin{bmatrix} N^s & 0 \\ N^{s'} & N^f \end{bmatrix} = T N, \qquad (2.16) $$

where $N^s$ and $N^f$ are stoichiometry matrices for the slow and fast part, respectively. The matrix $T$ has been chosen to match knowledge on the biochemical reactions involved. The variables $x^s \in \mathbb{R}_+^m$ typically represent total protein concentrations, obtained by summing the concentrations of free and complexed proteins. The variables $x^f \in \mathbb{R}_+^{n-m}$ include protein complexes and metabolites.

The introduction of the fast and slow variables leads to the following reformulation of the system (2.5):

$$ \frac{dx^s}{dt} = N^s v^s(x^s, x^f), \qquad x^s(0) = x_0^s, \qquad (2.17) $$

$$ \frac{dx^f}{dt} = N^{s'} v^s(x^s, x^f) + N^f v^f(x^s, x^f) \approx N^f v^f(x^s, x^f), \qquad x^f(0) = x_0^f \qquad (2.18) $$

where $v^s(x^s, x^f) \in \mathbb{R}^p$ and $v^f(x^s, x^f) \in \mathbb{R}^{q-p}$ are rate equations for the slow and fast system, respectively. The variable transformation ensures that slow variables change through slow reactions only. In system (2.18), however, we notice that the fast system still includes a slow time scale. This can be seen when we apply the linear transformation to our activation network model, which gives the system in Figure 2.1(c). The consumption rates of the three fast variables by growth dilution and protein degradation (represented by the parameter $\gamma$) are on a slow time scale compared to the rate of complex dissociation or export of cAMP. We neglect the slow term, that is, we assume that $N^{s'} v^s(x^s, x^f) \ll N^f v^f(x^s, x^f)$, as done in system (2.18). The principle was applied to the activation network model in Figure 2.1(b), which gives the slow-fast system in Figure 2.1(c).

Now that we have reformulated our system, we can apply the QSS assumption, based on the hypothesis that the fast variables instantaneously adapt to changes in the slow variables. This amounts to setting $dx^f/dt = 0$ in (2.18). From a mathematical point of view, the QSS approximation restricts the system dynamics to a manifold of lower dimension, an approximation of the so-called slow manifold (see [Heinrich and Schuster, 1996]). After an initial transient, the dynamics of the fast system can be well approximated by an algebraic function of the slow variables: $x^f = g(x^s), g: \mathbb{R}_+^m \to \mathbb{R}_+^{n-m}$. The following QSS model describes the dynamics of the system on the slow manifold:

$$ \frac{dx^s}{dt} = N^s v^s(x^s, g(x^s)), \qquad x^s(0) = x_0^s, \qquad (2.19) $$

The QSS approximation reduces the dimension of the system, but generally also the number of parameters, due to the fact that some of these can be lumped as they no longer independently
occur. The application of the QSS assumption to the model of the example network is shown in Fig. 2.1(d-e). We obtain an algebro-differential model with two state equations and one algebraic equation. The QSS approximation has reduced the total number of parameters from 19 to 13. Numerical simulation using physiologically relevant parameters shows how, after a transient time, the solution of the reduced model converges to the solution of the original one.

I have used this type of QSS approximation in various studies of the carbon starvation response network in E. coli [Baldazzi et al., 2010, 2012, Ropers et al., 2006, 2011], as well as transcription-translation [Belgacem et al., 2018], and mRNA degradation [Etienne et al., 2020] (Chapter 7.1). In [Ropers et al., 2011], for instance, we show by numerical studies that a QSS model of a large model of this system is a good approximation of the original model. The reduction allowed to simplify the original model with 14 variables and 63 parameters into a model with 9 variables and 59 parameters. The reduced model no longer includes parameters for which experimental data are generally not available, such as the association and dissociation rate constants, or model variables such as intermediary complexes. In Section 2.2.2, we will discuss how we can further simplify this type of model to qualitatively analyse its dynamics.

### 2.2.2 Model approximation by means of piecewise-linear functions

How to mathematically model biochemical networks when quantitative data is scarce? Even numerical simulation and analysis of models as simple as the one represented in Figure 2.1(e) can be complicated, due to a generalized lack of physiological values for parameters. Kinetic constants for gene expression like the maximal synthesis rate and the degree of cooperativity in the regulatory mechanisms are typically missing in most cases. Parameter estimation from experimental data could alleviate the problem, but this requires kinetic gene expression data that are also often missing. If no data can be acquired specifically for model estimation, an alternative is to further reduce the model and simplify its mathematical form.

During my postdoctoral research and in later collaborative projects, I focused on one specific type of simplification, the approximation of sigmoidal Hill functions by means of step functions [Glass and Kauffman, 1973, Mestl et al., 1995] (for review, see de Jong and Ropers 2006a,b). From a biological point of view, the use of step functions corresponds to the assumption that gene activity is switched on or off abruptly instead of progressively, when the concentration of the regulatory protein crosses a threshold. An example is shown in Figure 2.2 for two genes, *rrn* and *fis*, the transcription of which is regulated by protein Fis in a cooperative manner [see Ropers et al., 2006, and references therein]. We have modelled the rate of expression of stable RNAs in response to variations of the Fis concentration by a positive Hill function:

$$f_{rrn}(x_{fis}) = \kappa_{rrn} h^{+}(x_{fis}, \theta_{fis}^{1}, n_{1}), \text{ with } h^{+}(x_{fis}, \theta_{fis}^{1}, n_{1}) = \frac{x_{fis}^{n_{1}}}{x_{fis}^{n_{1}} + (\theta_{fis}^{1})^{n_{1}}} \quad (2.20)$$

where $\kappa_{rrn}$, $n_1$, and $\theta_{fis}^1$ are constants denoting the synthesis rate of stable RNAs, the Hill number and the dissociation constant of Fis, respectively. The function $f_{rrn}(x_{fis})$ implies that stable RNAs are maximally expressed at a rate $\kappa_{rrn}$ if $x_{fis} > \theta_{fis}^1$, whereas they are almost not expressed if $x_{fis} < \theta_{fis}^1$. The expression of $f_{rrn}(x_{fis})$ can be simplified by a step function, thus
<figure><img src="image_6.png" /><figcaption>**Figure 2.2** – Representation of the regulation of gene expression by means of step functions. Regulation of the expression of genes (a) *rrn* and (b) *fis* by the protein Fis. (c) Activity of the promoter *rrnP1* ($\alpha_{rrn}$), normalized on the scale 0 to 1, as a function of the concentration of transcriptional regulator Fis. (d) Idem for the activity of promoter *fisP* $\alpha_{fis}$. The dotted blue line represents the positive (c) or negative (d) step function used to approximate promoter activities modelled by Hill functions. Adapted from [Ropers et al., 2006].</figcaption></figure>

eliminating the often unknown cooperativity coefficient $n_1$ (Figure 2.2(c)):

$$f_{rrn}(x_{fis}) \approx \kappa_{rrn} s^{+}(x_{fis}, \theta_{fis}^{1}), \quad \text{with } s^{+}(x, \theta) \begin{cases} 0, & if \ x < \theta, \\ 1, & if \ x > \theta. \end{cases} \qquad (2.21)$$

In the same manner, the auto-inhibition of Fis synthesis can be modelled by a negative Hill function, which we subsequently approximate by a negative step function (Figure 2.2(d)):

$$f_{fis}(x_{fis}) = \kappa_{fis} h^{-}(x_{fis}, \theta_{fis}^{2}, n_{2}), \quad \text{with } h^{-}(x_{fis}, \theta_{fis}^{2}, n_{2}) = \frac{(\theta_{fis}^{2})^{n_{2}}}{x_{fis}^{n_{2}} + (\theta_{fis}^{2})^{n_{2}}}, \qquad (2.22)$$

and

$$f_{fis}(x_{fis}) \approx \kappa_{fis} s^{-}(x_{fis}, \theta_{fis}^{2}), \quad \text{with } s^{-}(x, \theta) \begin{cases} 1, & if \ x < \theta, \\ 0, & if \ x > \theta. \end{cases} \qquad (2.23)$$

The approximation can be extended to more complex systems, as in (2.19), where the dynamics of the fast system at steady state is determined by changes in the slow variables. The fast processes such as complex formation, metabolism, and signalling pathways introduce a coupling between slow variables, represented by a combinatorial expression of sigmoidal functions. For instance, following the application of the quasi-steady-state approximation in Figure 2.1(e), the term for the synthesis rate of protein CRP is a function of the CRP·cAMP concentration $f_{crp} = \kappa_c^1 + \kappa_c^2 h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m})$, which couples the CRP synthesis rate to the CRP and Cya concentrations because the steady-state concentration of the complex
$x_{c\sim m} = h^+(u_s, \theta_s, m_s) x_c x_y)/(K_4 K_3 + h^+(u_s, \theta_s, m_s) x_y)$ depends on these two variables. The resulting dependency of the CRP synthesis rate is shown normalized on the scale 0 to 1 in Figure 2.3(a).

Two steps are needed to approximate multivariate Hill functions like $h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m})$. In [Ropers et al., 2011], we showed that we can first approximate the multivariate function by a product of Hill functions:

$$h^+(g(x^s), \theta, n) \approx h^+(x_1^s, \theta_1, n_1) \times h^+(x_2^s, \theta_2, n_2) \times \dots \text{and:} \quad (2.24)$$

$$h^-(g(x^s), \theta, n) = 1 - h^+(x_1^s, \theta_1, n_1) \times h^+(x_2^s, \theta_2, n_2) \times \dots \quad (2.25)$$

We subsequently approximate the Hill functions by step functions:

$$h^+(g(x^s), \theta, n) \approx s^+(x_1^s, \theta_1) \times s^+(x_2^s, \theta_2, n_2) \times \dots \text{and:} \quad (2.26)$$

$$h^-(g(x^s), \theta, n) = 1 - s^+(x_1^s, \theta_1, n_1) \times s^+(x_2^s, \theta_2, n_2) \times \dots \quad (2.27)$$

Applied to the reduced model of the activation network in Figure 2.1(e), we obtain an approximation for the multivariate function $h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m}) \approx h^+(x_y, \theta_y, m_y) h^+(x_c, \theta_c, m_c) h^+(u_s, \theta_s, m_s)$, which we further approximate by a product of step functions $h^+(x_{c\sim m}, \theta_{c\sim m}^1, m_{c\sim m}) \approx s^+(x_c, \theta_c^1) s^+(x_y, \theta_y^1) s^+(u_s, \theta_s)$. This eventually gives the so-called piecewise-linear (PL) model in Figure 2.3(a). The model describes the regulatory logic of CRP expression: the protein is expressed if glucose is depleted ($u_s > \theta_s$), and if there is enough enzyme Cya to produce cAMP ($x_y > \theta_y$) and transcription factor CRP ($x_c > \theta_c$). We showed by numerical studies in [Ropers et al., 2011] that PL models are good approximations of reduced and original models. In the case of the carbon starvation network studied in this paper, the PL approximation reduces the QSS model of 59 parameters to a PL model of 46 parameters [Ropers et al., 2011].

<figure><img src="image_7.png" /><figcaption>**Figure 2.3** – (a) Approximation of complex sigmoidal functions by a product of Hill functions. The surface represents the plot of $h^+(x_{c\sim m}, \theta_{c\sim m}, m_{c\sim m})$ as a function of $x_y$ and $x_c$. It is fitted with the product of Hill functions $h^+(x_y, \theta_y, m_y) h^+(x_c, \theta_c, m_c) h^+(u_s, \theta_s, m_s)$. $\theta_c$ and $\theta_y$ denote the threshold values determined by curve fitting. (b) PL model for the activation network. From [Ropers et al., 2011].</figcaption></figure>

### 2.2.3 Discussion and perspectives

In this chapter, I presented work related to the approximation of nonlinear ODE systems, allowing to develop a large variety of biochemical network models, from more detailed to more
abstract models. The various reduced models are rigorously derived, which results in good approximations that preserve the main features of the global system dynamics [Ropers et al., 2011]. Quasi-steady-state and piecewise-linear approximations will be repeatedly used in the various biological problems addressed in the following chapters.

A key step of the method presented in Section 2.2.1 is the ability to solve a fast system at steady state. This is no longer possible when biochemical networks are much larger and include slow variables that are indirectly coupled through numerous intermediate metabolic and signalling processes. This point was addressed in follow-up studies in the context of the post-doctorate of Valentina Baldazzi, with the development and application of a method allowing the inference of PL models from more complex QSS models by means of time-scale arguments and sensitivity criteria from metabolic control analysis [Baldazzi et al., 2010, 2012]. I also worked on time-scale separation in the context of transcription-translation models in collaboration with the team of Jean-Luc Gouzé (Inria Sophia Antipolis - Méditerranée) [Belgacem et al., 2018].

Other related work was carried out in the context of the PhD thesis of Stefano Casagranda, whom I co-supervised with Jean-Luc Gouzé. We developed an approach, Principal Process Analysis, to determine the contribution of each biological process to the output of a dynamical system. Due to the general form of biochemical network models in (2.5), the rates of the processes appear in a linear additive manner in the right-hand side of the ODEs. In [Casagranda et al., 2015, 2018], we introduced new quantities to weigh the influence of each process on the dynamical change of each model variable. Processes that are inactive in some time windows, because they do not influence the dynamics of the system, are removed from the model. This allows the creation of submodels for each time window that only contain the active processes. This procedure leads to the simplification of the system to its core mechanisms. The simplified system can be further studied, to understand the role of each active process in the system dynamics. Formally, the method is not a model reduction approach because it does not preserve mass balances. However, it allows to dissect the complex dynamics of original models through the analysis of simplified versions of these models in given time windows.

There is no firm response to the question as to which model formalism is the most appropriate. A trade-off has to be found between the quantitative and mechanistic information available and the complexity of the problem studied. ODE models of gene expression and metabolism and their reduced QSS models offer a powerful framework if molecular concentrations and kinetic constants are known or can be precisely inferred from experimental data. When these models are too large for a dynamical analysis, or when information on gene expression dynamics, signalling or metabolic pathways is incomplete or missing, PL models are a good alternative. This is the situation faced in the following chapter with the study of the carbon starvation response of the bacterium *E. coli*. The PL model of the gene regulatory network controlling the stress response is analysed by means of qualitative approaches.
<figure><img src="image_8.png" /></figure>
# Qualitative analysis of the dynamics of gene regulatory networks

The work described in this chapter represents the main achievement of my post-doctorate in the Helix Bioinformatics group [Ropers et al., 2006]. The piecewise-linear (PL) model developed in this study is the first attempt to connect *E. coli* growth with a network of transcription regulators involved in the adaptation of *E. coli* to its nutritional environment. At the time, there was a lot of information in the literature about these regulators, but expression data was sparse, often obtained in diverse steady-state environmental conditions, and for different strains. Because of the lack of quantitative information, I used a qualitative modelling and simulation approach to analyse the functioning of the network of transcription regulators.

In the following section, I will illustrate this approach with the PL model of a cross-inhibition network. This simple network is part of the carbon starvation network, whose dynamics will be analysed in Section 3.2.

## 3.1 Qualitative modelling and simulation of piecewise-linear models

An advantage of the use of step functions is that they facilitate the analysis of the qualitative dynamics of the PL models (de Jong et al. 2004; for review, de Jong and Ropers 2006a,b). The threshold values of step function variables partition the phase space into hyper-rectangular regions, in each of which the system behaviour is qualitatively homogeneous. The continuous phase-space dynamics of the system can be discretised into a state transition graph. This graph is composed of states corresponding to regions of the phase space as well as transitions between these states. The state transition graph describes the possible qualitative behaviours of the system and allows to determine attractors and their attainability. Such analysis is difficult to carry out with the corresponding quantitative QSS models [Ropers et al., 2007].

An example of qualitative analysis is shown with a small cross-inhibition network in Figure 3.1(a). The gene *crp*, which we have encountered in previous examples, is involved in a positive feedback loop with gene *fis*, while *fis* expression is also limited by autoinhibition. Panel (b) of the figure shows the corresponding PL differential equation model, in the case of carbon starvation (when CRP has already been activated by the binding of cAMP upon glucose depletion). In the absence of quantitative information, we qualitatively order the concentrations and threshold parameters (panel (c)). For instance, the concentration of protein Fis varies between a minimal and a maximal value denoted 0 and $max_{fis}$, respectively. When its concentration reaches the threshold $\theta_{fis}^1$, the protein binds to *crp* promoter and inhibits CRP expression. At
higher intracellular levels, when it reaches the threshold concentration $\theta_{fis}^2$, the protein represses its own expression by binding to its promoter region. The dynamics of the system in the phase plane is shown in panel (d). The system possesses three equilibrium points, two stable points – one characterized by a high concentration of Fis and low concentration of CRP, and the other by a low concentration of Fis and high concentration of CRP – and one unstable point. In the upper left region characterized by $0 < x_{fis} < \theta_{fis}^1$ and $\theta_{crp} < x_{crp} < \max_{crp}$, the step function $s^-(x_{crp}, \theta_{crp})$ evaluates to 0, while $s^-(x_{fis}, \theta_{fis}^1)$ and $s^+(x_{fis}, \theta_{fis}^2)$ are equal to 0. The PL model in this specific region simplifies to the following system:

$$ \frac{\dot{x}_{fis}}{dt} = -\gamma_{fis} x_{fis}, \qquad (3.1) $$

$$ \frac{\dot{x}_{crp}}{dt} = \kappa_{crp} - \gamma_{crp} x_{crp}. \qquad (3.2) $$

Trajectories in this region asymptotically converge to the point $x_{fis} = 0$ and $x_{crp} = \frac{\kappa_{crp}}{\gamma_{crp}}$. To each region of the phase space corresponds a simplified system of linear equations with an associated dynamical behaviour. After discretization of the phase space, we obtain the state transition graph in panel (d), which includes three qualitative equilibrium states, two stable and one unstable one. As can be seen with this simple example, the qualitative modelling approach allows a quick scan of the qualitative dynamics of the system, without numerical information on parameter values [de Jong et al., 2004]. It has been implemented in the computer tool Genetic Network Analyzer [de Jong et al., 2003, Batt et al., 2012], allowing the analysis of large state transition graphs of complicated gene expression network models.

## 3.2 Qualitative analysis of the carbon starvation response in E. coli

A variety of processes are involved in the adaptation of E. coli bacteria to a carbon source, from the remodelling of gene expression, metabolism, and DNA topology to the adaptation of the growth rate, and the lack of mechanistic information. Based on literature data, we reconstructed a first network of six genes that are believed to play a key role in the response of the cell to carbon source availability (Figure 3.2(a)). The network includes genes involved in the transduction of the carbon starvation signal (the global regulator of transcription crp and the adenylate cyclase cya; panels (a) and (b)), metabolism (the global regulator fis), cellular growth (the rrn genes coding for stable RNAs, needed in high numbers in exponential phase of growth for ribosome biogenesis and protein synthesis), and DNA supercoiling, an important modulator of gene expression (the topoisomerase topA and the gyrase gyrAB; panels (a) and (c)). The resulting network is of course far from exhaustive and leaves aside other known regulators. However, it allows to analyse how the key genes included function together. The outcome of their interactions is hard to predict because they are involved in many feedback loops, such as the mutual inhibition of fis and crp, the auto-inhibition of fis mediated by Fis itself and by the DNA supercoiling or the auto-inhibition of CRP and Cya. Do these network components and shared interactions allow to reproduce and explain biological observations?

The network was reconstructed through a bottom-up approach, by gathering literature and database knowledge on the different genes. As described in Chapter 2, writing the correspond-
<figure><img src="image_9.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_10.png" /></figure>

$$
\begin{align*}
\dot{x}_{\text{fis}} &= \kappa_{\text{fis}} s^{-}(x_{\text{crp}}, \theta_{\text{crp}}) s^{-}(x_{\text{fis}}, \theta_{\text{fis}}^{2}) - \gamma_{\text{fis}} x_{\text{fis}}, \\
\dot{x}_{\text{crp}} &= \kappa_{\text{crp}} s^{-}(x_{\text{fis}}, \theta_{\text{fis}}^{1}) - \gamma_{\text{crp}} x_{\text{crp}},
\end{align*}
$$

(b)

$$
\begin{align*}
0 < \theta_{\text{fis}}^{1} < \theta_{\text{fis}}^{2} < \max_{\text{fis}}, \quad \theta_{\text{fis}}^{2} < \kappa_{\text{fis}}/\gamma_{\text{fis}} < \max_{\text{fis}}, \\
0 < \theta_{\text{crp}} < \max_{\text{crp}}, \quad \theta_{\text{crp}} < \kappa_{\text{crp}}/\gamma_{\text{crp}} < \max_{\text{crp}}.
\end{align*}
$$

(c)

<figure><img src="image_11.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e)</figcaption></figure>

**Figure 3.1** – (a) Simple genetic regulatory network composed of the genes *fis* and *crp*. In conditions of carbon starvation, when CRP is activated by cAMP, Fis and CRP inhibit each other’s expression. (b) PL differential equation model of the small network. (c) Inequality constraints on threshold and rate parameters associated to PL models. (d) Sketch of the dynamics in each domain of the phase space for the two-gene network. Dots represent the equilibrium points of the system. (e) State transition graph for the two-gene network. QS denotes a qualitative state. The qualitative equilibrium states are circled. Adapted from [Ropers et al., 2007].
<figure><img src="image_13.png" /><figcaption>**Figure 3.2** – (a) Network of key genes, proteins, and regulatory interactions involved in the carbon starvation response network in *Escherichia coli*. The boxes ‘CRP activation’ and ‘supercoiling’ are detailed in Panels (b) and (c). The graphical conventions [Kohn, 2001] are explained in the legend. (d) PL differential equation and parameter inequality constraints for the topoisomerase TopA. From [Ropers et al., 2006].</figcaption></figure>
ing PL model required the prior development of a kinetic model and its reduction based on time-scale separation to obtain a QSS model describing the coupling of slow variables – Fis, CRP, GyrAB, TopA, stable RNAs – with the fast processes – DNA supercoiling and activation of CRP by the phosphotransferase system. The coupling is schematically represented in Figure 3.2(b,c). We eventually obtained a PL model of seven variables, one concentration variable for each of the six proteins and one input variable representing the presence or absence of the carbon starvation signal. Seven differential equations, one for each variable, and forty inequality constraints describe the dynamics of the system. As an illustration, the differential equation and the parameter inequality constraints for the state variable $x_{topA}$ are given in Figure 3.2(d). For instance, the constraint $0 < \frac{\kappa_{topA}^1}{\gamma_{topA}} < \theta_{topA}^1$ expresses that without stimulation of the topA promoter, the TopA concentration decreases towards a background level, below the threshold $\theta_{topA}^1$.

The PL model is able to predict growth arrest and the entry into stationary phase following carbon depletion, as well as growth resumption to exponential growth when glucose is available again [Ropers et al., 2006]. Two state transition graphs are obtained for each of these simulation scenarios. One includes all possible qualitative behaviours from the initial state representative of the exponential growth conditions to the equilibrium state corresponding to stationary phase conditions and the opposite for the other graph. As an illustration, a representative path along the graph for carbon starvation response is shown in Figure 3.3. It reproduces experimental observations from the literature, such as the decrease of protein Fis concentration in response to carbon depletion and growth arrest due to the reduced production of stable RNAs. The chain of molecular events can be deduced from the prediction and the network structure. The lack of glucose activating CRP, the protein represses Fis synthesis. The decrease of the protein concentration alleviates the inhibition of CRP synthesis. The signal of carbon starvation is thus amplified by the mutual inhibition module involving *fis* and *crp*. The reduction of Fis levels affects the DNA supercoiling module as well as cell growth is halted by the arrest of stable RNA production. Predictions of growth resumption when glucose is added were more surprising, as damped oscillations of Fis, stable RNA, and GyrAB concentrations were predicted to occur [Ropers et al., 2006]. This led to launching experiments in the group to verify the reality of this behaviour. I will come back to this in Chapter 5.

## 3.3 Discussion and perspectives

The PL model developed in this study is the first attempt to connect *E. coli* growth with a network of transcriptional regulators based on the literature information. The qualitative analysis of the PL model allowed to bypass the lack of quantitative information and understand the cascades of molecular processes involved in the stress response [Ropers et al., 2006]. This work also allowed to point at inconsistencies in the literature concerning the regulation of transcription factors. We will come back to this point in Chapter 5, where the global control of gene expression will be shown to be an important regulatory mechanism of the environmental adaptation of *E. coli*.

Although PL models are reduced and approximated versions of larger nonlinear ODE models, they can be complex and generate large state-transition graphs including many
<figure><img src="image_14.png" /><figcaption>**Figure 3.3** – Temporal evolution of the protein and stable RNA concentrations in a typical qualitative behaviour in the state transition graph, when Signal is present ($\theta_s < u_s \le \max_s$). The behaviour represents the molecular events accompanying the transition from exponential to stationary phase following carbon starvation. From [Ropers et al., 2006].</figcaption></figure>
alternative qualitative behaviours. This makes it impossible to analyse the system dynamics by manual inspection and to compare model predictions to available experimental data. Ideally, one would like to query if the graph is consistent with our understanding or with the partial and heterogeneous information from the literature. For instance, is there a path in the graph reproducing the experimental observation that Fis concentration decreases when there is a lack of glucose? And is there a concomitant increase of the CRP concentration? This type of problem was addressed in subsequent work by Grégory Batt and Pedro Monteiro during their PhD theses in our group. The qualitative modelling approach was reformulated to allow the application of formal verification tools based on model checking, together with the development of an approach for the automatic translation of frequent biological questions into the temporal logic formulas used in model checking. This facilitates the application of model-checking approaches for the analysis of state transition graphs by non-expert users. I participated in the development of these queries and the application of model-checking to analyse the PL model of the network in Figure 3.2(a) [Batt et al., 2005] and an extension of this network that I derived [Monteiro et al., 2008].

The models that we developed were reconstructed and manually parametrised from the literature information. I collaborated on several other projects in this manner, for the piecewise-linear modelling of the transcriptional regulatory network controlling the response of the model eukaryote *Saccharomyces cerevisiae* to the agricultural fungicide mancozeb [Monteiro et al., 2011], or the network involved in the adaptation of food pathogens to osmotic stress [Ropers and Métris, 2016, Métris et al., 2017]. Alternative approaches exist though, such as model inference from reporter gene data as done in our group in the framework of the PhD thesis of Diana Stefan [Stefan et al., 2015]. Note that the PL models can be further approximated by logical models [Abou-Jaoudé et al., 2016, de Jong, 2002, Kauffman, 1969, Thomas and d'Ari, 1990]. This opens the door to other advanced tools and methods for the analysis of discrete models. In [Corblin et al., 2009], for instance, I collaborated with the group of Laurent Trilling (Univ. Grenoble Alpes) on the reconstruction of discrete versions of the model in [Ropers et al., 2006] that are consistent with biological data expressed as constraints.

One may argue that not everyone masters model reduction into QSS models and their approximation into PL models whose dynamics can be qualitatively analysed. Based on my experience, I advised the (former) company Genostar in a follow-up project. The company was involved in the commercialization and distribution of a bioinformatics platform including Genetic Network Analyzer. The project led to the development of an approach facilitating the rigorous building of PL models able to address the problem of coupling slow variables through fast processes, which was implemented as a graphical tool in subsequent releases of Genetic Network Analyzer [Batt et al., 2012]. This approach facilitates PL model building by less informed modellers.

The work described in this chapter was motivated by the lack of quantitative information to parametrize the PL model and verify its predictions. This led Hans Geiselmann and his group to mount an experimental program, in which I was involved as well, for the development of strains and plasmids expressing reporter genes that could be used for monitoring gene expression in
various strains and conditions. This work was in line with the general trend of biology becoming increasingly quantitative. As a consequence of these developments, there was an interest in trying alternative approaches based on the analysis of quantitative ODE models or their reduced versions. I will describe this work in Chapters 6 to 7. The availability of new experimental approaches with new types of data opens new methodological questions as well, such as the robust inference of biological quantities from the primary data than can be used for the analysis and validation of quantitative ODE models. These methodological questions will be the subject of the following chapter.
# Analysis of dynamical gene expression and metabolomics data

In the course of the years, biology has moved from experimental approaches based on qualitative or semi-quantitative experiments (e.g., western-blots) to more quantitative ones. The most recent ones include transcriptomics and proteomics data for instance, measurements of gene expression along growth by means of reporter genes, or dynamical metabolomics measurements. The two later approaches prove essential in a number of modelling studies described in the following chapters. Most of the time, the time-series data need to be processed in order to estimate the relevant biological quantities. For purposes such as model validation against experimental data, parameter estimation, or model-based analysis of the data, it is important that these estimations be unbiased and robust to experimental noise, while remaining able to capture precisely rapid variations of the signals, which commonly occur during growth-phase transitions.

In this chapter, I will present work on the use of regularized linear inversion methods for the reconstruction of interesting quantities from metabolomics (Section 4.1) or gene expression (Section 4.2) data. I chose to highlight these papers because of their relevance for many laboratories. Indeed time-series data accumulate in these places because they are becoming easier to acquire, but their information content is far from exploited due to the difficulties to analyse them. The method for gene expression analysis was developed by Valentin Zulkower as part of his PhD thesis, which I co-advised with Hidde de Jong and Hans Geiselmann [Zulkower et al., 2015]. The method for metabolomics data analysis was developed by Eugenio Cinquemani. For this study, I brought the problem, analysed the data, and coordinated the project [Cinquemani et al., 2017]. The method and results were presented at the joint 2017 ISMB-ECCB conference, and published in the corresponding special issue of *Bioinformatics* [Cinquemani et al., 2017].

## 4.1 Estimation of time-varying growth, uptake and secretion rates from dynamic metabolomics data

The accumulation of extracellular metabolites or their disappearance from the growth medium provides interesting information about intracellular physiology [Kell et al., 2005]. The time profiles are used to compute uptake and secretion rates that can be related to intracellular fluxes in flux balance and metabolic flux analyses, as will be shown in Chapter 6. Despite its apparent simplicity, the problem of estimating time-varying uptake and secretion rates from measurements of extracellular metabolites is challenging. First the available data are noisy and sparse, despite continuous progress in metabolomics methods. Second, the time-course
<figure><img src="image_15.png" /><figcaption>**Figure 4.1** – Reconstruction of growth and metabolite exchange rates from E.coli diauxic growth experiments with the method in [Cinquemani et al., 2017]. (a) Model scheme defining variables and rates. (b) Data (circles) from [Morin et al., 2016] and EKS estimates with 95% credibility intervals (solid curve and shaded band) of biomass and concentration profiles at all times. (c) EKS rate estimates and credibility intervals (red curve and shaded band). Vertical blue lines show the periods of fast transitions detected in data-preprocessing).</figcaption></figure>

profiles of different extracellular metabolites are strongly correlated. Uptake and secretion rates are indeed proportional to the size of the growing population of cells consuming or producing the metabolites [Stephanopoulos et al., 1998]. Third, the profiles are subject to discontinuities, resulting from sudden changes in the functioning of metabolism. Transition from growth on glucose to acetate is one example, see Figure 4.1(b) for an example (data from Chapter 6; [Morin et al., 2016]). The acetate suddenly stops accumulating in the growth medium upon glucose depletion and is right after consumed by cells as they resume their growth. Such sharp transition of metabolic regime is difficult to account for. The consequence is that exometabolomics data are often under-exploited, in the sense that in many studies the time-varying rates of uptake and secretion are not computed to only focus on the steady-state regime for which data analysis is easier. More advanced approaches are needed to analyse the data.

ODE models of microbial growth are traditionally used to estimate exchange/growth rates from metabolomics and growth data. The general form of the model describes the dynamics of biomass concentration $b(t)$ and concentration of metabolite $i$ at time $t$ $c_i(t)$ [de Jong et al., 2017a]:

$$ \begin{aligned} \dot{b}(t) &= \mu(t) b(t), \\ \dot{c}_i(t) &= r_i(t) b(t), \quad i = 1, \dots, n. \end{aligned} \qquad (4.1) $$

The model is based on the assumption that changes in concentrations $c_i(t)$ are due to the growth, uptake and secretion rates, omitting the negligible degradation rate of extracellular metabolites, inflow and outflow of the medium in the bioreactor. Growth rate $\mu(t)$ and exchange
rates $r_i(t)$ must be reconstructed over a whole experimental period $t \in \mathcal{T}$ from noisy measure-
ments of $b$ and of $c_i$ taken at possibly different time instants. Traditionally, this is achieved by
differentiating a fit of the corresponding measurements and computing the estimated rates by
means of Equation 4.2:

$$
\begin{align}
\mu(t) &= \frac{\dot{b}(t)}{b(t)}, \tag{4.2} \\
r_i(t) &= \frac{\dot{c}_i(t)}{b(t)}. \nonumber
\end{align}
$$

Such approach is said indirect, as the data are first smoothed before the quantities of interest
are reconstructed via the model in Equation 4.2. This results in the propagation of estimation
errors that cannot be controlled. While spline fitting of the data is often the preferred smoothing
method [Wahba, 1990], the optimal placement of spline knots is a difficult problem, in particular
when there are sudden changes of regime, and the quality of the estimate needs to be assessed
in further processing, typically through bootstrap approaches. Finally, the approach does not
account for the coupling of the different model equations through the biomass, while the ex-
ploitation of this important information could improve the rate estimation procedure. Instead
each rate is estimated independently from the others.

We addressed the above issues in [Cinquemani et al., 2017]. Our (direct) approach is
based on the use of the explicit model in Equation 4.1, together with an appropriate statistical
approach for the estimation of the exchange rates. We considered the reconstruction of rates
from measured concentrations and biomass as an input estimation problem, which can be
solved using linear inversion methods. That is, the growth, uptake and secretion rates are
the input (u(t) = [μ(t) r₁(t) ... rₙ(t)]ᵀ), and the biomass and metabolite concentrations the
output (x(t) = [b(t) c₁(t) ... cₙ(t)]ᵀ) of the linear system dx(t)/dt = u(t)x(t), so that rates
can be estimated by linear inversion from the noisy biomass and metabolite concentration
measurements y(t) modelled as: y(t) = x(t) + e(t), t ∈ ℐ, with e(t) a Gaussian measurement
noise. The input profile u(t) is assumed to be piecewise continuous, so that the solution of the
ODE system is well determined, but not necessarily smooth.

As it stands, the problem of estimating uptake/secretion and growth rates is ill posed. There
are an infinity of solutions that explain the data for a given set of initial conditions, some of which
are not realistic from a biological point of view. For instance, irregular ("wiggly") profiles may
fit slowly changing measurements of the biomass or the metabolite concentrations. The problem
must be regularized to obtain a unique, acceptable solution [Wahba, 1990]. While it avoids over-
fitting the data, regularization is challenging because the data are characterized by a combination
of slow and fast variations of metabolite concentrations and biomass. In [Cinquemani et al.,
2017], we solved this by formulating a Bayesian regularized estimation problem. Each unknown
rate profile uᵢ(t) with i = 1, ..., n + 1 is modelled as the outcome of a random Gaussian process
dvᵢ(t) / dt = γ(t) ωᵢ(t) and duᵢ(t) / dt = vᵢ(t), where wᵢ is the standard white Gaussian noise and γᵢ(t) the
regularization factor. uᵢ(t) is hence modelled as a double-integral of white noise, which implies it
being continuously differentiable with variability (the probability distribution of its derivative)
determined by the magnitude of γᵢ(t) > 0. γᵢ(t) is a function of time, where larger values of the
factor around a time point allow for rapid changes of uᵢ(t) around that time. By this modelling,
the problem of estimating uᵢ at any time t becomes that of computing the best estimates in a
Bayesian perspective, provided the measurement model for the set of all data ℳ = ℳ₁, ..., ℳₙ₊₁
with $\mathcal{Y}_i = y_i(t) : t \in \mathcal{T}_i$, to calculate the a posteriori expectation

$$ \mathbb{E}[(r_1(t), \dots, r_n(t), \mu(t))|\mathcal{Y}] \qquad (4.3) $$

at all times of interest $t \in \mathcal{T}$. This approach favours smooth solutions. Two data preprocessing steps are carried out separately to obtain Bayesian priors for the smoothing factors for the slow and fast dynamics, and the detection of regions where fast dynamics take place. The first one is based on a cubic spline interpolation of the data and generalized cross validation, separately for every $i$, while the second detects times at which concentrations drop to zero, a fact generally associated with a change of metabolic regime. The regularization estimation problem is solved by a dynamical smoothing approach based on an Extended Kalman Smoother (EKS) [Kailath et al., 2000]. It allows to obtain a smoothed estimate for both the rates and states (concentrations and biomass) at time $t$, along with credible intervals, given past and future measurements. In practice, the Kalman filter - in its extended version adapted for non-linear systems - allows in a first step to estimate states and rates at time $t$ from the knowledge of past measurements from 0 to $t$ and updates the solution given the measurement at $t+1$. In a second step, the smoothing allows to estimate states and rates on the interval $[0, T]$, given measurements $y(0: T)$. Smoothed estimates are more accurate than the filtered ones because more data are used [Cinquemani et al., 2017].

We tested the approach on simulated data and real data: those obtained from the batch cultures of wild-type *E. coli* obtained during diauxic growth on glucose and acetate (Chapter 6; [Morin et al., 2016]), as well as fed-batch culture data of *L. lactis* obtained at TBI [Cinquemani et al., 2017]. Results obtained with *E. coli* are shown in Figure 4.1. The method is able to capture the abrupt regime changes, while providing us with stable growth and exchange rates during periods of slow dynamics (either during growth on glucose before time 0 of glucose exhaustion or during growth on acetate, after time 0). We validated our approach by comparing metabolic flux analysis results obtained either with exchange/growth rates obtained from smoothing splines (those used in the data preprocessing step) or the EKS approach. The smoothing spline estimates yield distributions of intracellular fluxes that are much less precise and non-intuitive. On the contrary, the EKS approach allows to make accurate predictions of intracellular fluxes when compared to literature data. This shows the importance of precise rate estimation for metabolic flux analysis using the EKS approach. I used the approach for this precise reason in further studies of *E. coli* metabolism, such as the one on the role of glycogen in metabolic adaptation of *E. coli* described in Section 6.2.

## 4.2 Estimation of promoter activities and protein concentration profiles from reporter gene data

Reporter gene experiments, which enable to monitor gene expression at high time resolution in a non-intrusive way, are commonly used for the study of gene activities in response to environmental changes, such as in Chapter 5. In these experiments, a gene coding for a fluorescent (or luminescent protein) is designed to have the same promoter as a natural gene of the bacteria, and therefore be driven by the same regulations. A bacterial strain carrying this synthetic gene, either on the chromosome or on plasmid, is grown in a microplate or a
<figure><img src="image_16.png" /></figure>

**Figure 4.2** – Comparison of metabolic flux analysis results using rate estimates obtained with the EKS method in [Cinquemani et al., 2017] or smoothing splines. (a) Intervals of fluxes in the optimal solutions for 15 selected reactions in the carbon central metabolism obtained with the smoothing spline (grey) and EKS (red) estimates. (b) Comparison of predicted and measured intracellular fluxes for eleven reactions in the carbon central metabolism (EX<sub>glc</sub>, PGI, G6PDH2r, PGMT, GAPD, PYK, PDH, PPC, MDH, TKT1, and TKT2). Fluxes are given relatively to the specific glucose consumption rate.
bioreactor with frequent measurements of fluorescence (or luminescence) and absorbance, which are then processed using models of gene expression to estimate the temporal profiles of the population's growth rate, promoter activity (transcription rate) and intracellular concentration of the protein encoded by the gene of interest. The relation between promoter activity and the observed fluorescence and absorbance signals is indirect. Models of gene expression are thus needed to infer the quantities of interest. In past studies, we addressed the problem of estimating promoter activities through an indirect approach, where we first smoothed the primary data using cubic splines and then injected the resulting approximating functions into a gene expression model (see [de Jong et al., 2010] for the method and [Boyer et al., 2010] for its software implementation). As pointed out in Section 4.1, however, this leads to estimation errors that are difficult to control. In addition, and in the same way that sharp metabolic shifts are observed during growth transitions, gene expression abruptly changes as well, which is difficult to capture with smoothing splines. Here also, the problem of reconstructing interesting quantities from reporter gene signals requires advanced methods. In a subsequent study [Zulkower et al., 2015], we showed how the analysis of reporter gene data could be also formulated as an input estimation problem and be solved by means of regularized linear inversion.

I briefly illustrate the principle of the method with a simple ODE model describing reporter gene expression in one step:

$$ \frac{dR(t)}{dt} = a(t) V(t) - \gamma_r R(t), \qquad (4.4) $$

where $R(t)$ [mmol] is the time-varying quantity of the reporter protein in the growing bacterial population and $a(t)$ [mmol $min^{-1}$ $L^{-1}$] the synthesis rate of the reporter protein per unit of population volume $V(t)$ [L]. The latter is also called gene activity or promoter activity in the literature [Ronen et al., 2002], assuming a proportionality between transcription and translation rates. The fact that we consider total amount of molecules instead of concentrations as done usually allows to drop the term of growth dilution from the equation. $\gamma_r$ is the degradation constant of the reporter [$min^{-1}$] and can be easily measured [de Jong et al., 2010]. More detailed models can be used, for instance to describe the synthesis of the reporter mRNA, as well as the immature and mature forms of the reporter protein [Zulkower et al., 2015]. This is appropriate in situations where the fluorescent protein has a relatively long maturation time or the reporter mRNA is stable, for instance.

Absorbance and fluorescence measurements carried out in microplate readers are generally assumed to be proportional to the volume of the growing bacterial population $V$ and the total amount $R$ of reporter protein in the population, respectively. The following measurement model thus relates the absorbance $\tilde{V}$ and fluorescence $\tilde{R}$ measurements at time $t_i$ to the volume and reporter protein quantities:

$$ \begin{aligned} \tilde{V}(t_i) &= \alpha V(t_i) + \nu_i, \\ \tilde{R}(t_i) &= \beta R(t_i) + \nu'_i, \end{aligned} \qquad (4.5) $$

where $\nu_i, \nu'_i$ represent the measurement noise and $\alpha, \beta$ are unknown proportionality coefficients. The concentration of the reporter protein can be simply obtained from the ratio $\tilde{R}(t_i)/\tilde{V}(t_i)$.

We solve a first estimation problem, by inferring the growth rate of the population from the
<figure><img src="image_17.png" /><figcaption>**Figure 4.3** – Fluorescence and absorbance data obtained from reporter gene experiments in *E. coli* and estimations of growth rate, promoter activity and protein concentration from these data. The measured fluorescence and absorbance signals are shown in the top row. The estimations of growth rate, promoter activity and protein concentration are denoted by $\hat{\mu}(t)$, $\hat{a}(t)$, and $\hat{p}(t)$, respectively. The fluorescence signal, $\hat{a}(t)$, and $\hat{p}(t)$ have been divided by their mean as they have different orders of magnitude for *E. coli* genes *fis*, *gyrA*, *crp* and *acs*. For each signal, four replicates are shown, corresponding to different wells of the microplate.</figcaption></figure>

absorbance measurements via the following growth equation:

$$ \frac{d(\alpha V)(t)}{dt} = \alpha V(t) \mu(t) \approx \bar{V}(t) \mu(t), \qquad (4.6) $$

where $\mu(t)$ [min⁻¹] represents the growth rate and $\tilde{V}(t)$ an interpolated version of the measurements $\tilde{V}(t)$. As it is posed, the estimation of growth rate $\mu(t)$ and initial conditions $(\alpha V)(t_0)$ is a linear estimation problem, where $\mu(t)$ is the input (assumed to be piecewise constant) and $\tilde{V}(t)$ the observed output. Similarly to the analysis of metabolomics data, the problem is ill-posed and regularization is needed to solve the linear inversion problem. In the present case, we imposed a Tikhonov regularization on the first derivative of the growth rate [Zulkower et al., 2015], penalizing rapid successive variations, and we set the regularization parameter by generalized cross validation. We illustrate the results obtained in Figure 4.3 with data obtained from *E. coli* cells growing on glucose. Cells were transformed with plasmids carrying a transcriptional fusion of *gfp* with promoters of genes involved in bacterial adaptation to carbon source availability. The estimated growth rates are shown on the second row of Figure 4.3 for four replicates. As can be seen, the growth rate drops from its maximal value to zero upon glucose exhaustion. The estimation is thus able to capture the sudden arrest of bacterial growth.

Along with an estimated growth rate, the linear inversion gives us the estimated volume of the growing cell population $\tilde{\alpha}V(t)$, which we use in replacement of the cell volume in Equation 4.4
to solve the second estimation problem:

$$ \frac{dR(t)}{dt} = \frac{a(t)}{\alpha} \widehat{V(t)} - \gamma_r R(t). \qquad (4.7) $$

This is again a linear, time-varying system with input $a(t)$ and output $R(t)$, which can be solved by means of the same linear inversion approach as above. The reconstructed promoter activities for the four reporter genes in E. coli are shown in the third row of Figure 4.3. The method correctly infers the known fast changes in gene expression from the data, such as induction of *fis* following a nutrient upshift [Ali Azam et al., 1999] and *acs* upon glucose exhaustion [Berthoumieux et al., 2013, Wolfe, 2005], while avoiding over-fitting outside the transition region.

The host protein and its reporter counterpart share the same promoter and control sequences, but do not necessarily have the same degradation rate. Based on the knowledge of the degradation rate constant, it is also possible to reconstruct the concentration of the protein of interest expressed from the natural gene on the chromosome. To that aim, we develop the following model for the expression of the host protein:

$$ \frac{dP(t)}{dt} = \frac{a(t)}{\alpha} \widehat{V(t)} - \gamma_p P(t) \qquad (4.8) \\ p(t) = \alpha \frac{\widehat{P(t)}}{\widehat{V(t)}} $$

where $P(t)$ [mmol] is the time-varying quantity of the host protein, $p(t)$ [mmol L⁻¹] its concentration, and $\gamma_p$ [h⁻¹] its degradation constant. We assume $\gamma_p$ to be approximately known, bearing in mind that most bacterial proteins are stable, with half-lives well over 10 h [Larrabee et al., 1980]. If the reporter protein is stable as well, the default choice of $\gamma_p = \gamma_r$ usually leads to good results, in the sense of returning a protein concentration that is a smoothed version of the reporter concentration. We have again a linear system in Equations 4.7-4.8, where the equations for the reporter and host protein quantity are coupled by a shared input $a(t)$. The linear inversion procedure was adapted for linear equations with a shared input in [Zulkower et al., 2015], resulting in an estimate of $p(t)$ from absorbance and fluorescence measurements. The estimated protein concentrations for the E. coli reporter genes are shown in the bottom row of Figure 4.3. The degradation constant of Fis was measured in a former study ($d_p = 0.0065 \text{ min}^{-1}$; [de Jong et al., 2010]), whereas the other proteins were assumed long-lived ($d_p = 0.001 \text{ min}^{-1}$) like most E. coli proteins [Larrabee et al., 1980]. For instance, the estimated profile of Fis reproduces correctly the known transient accumulation of the protein following a nutrient upshift [Ali Azam et al., 1999, de Jong et al., 2010], which is consistent with its key role in activating the synthesis of stable RNAs needed for growth [Dennis et al., 2004].

More generally we tested the overall approach on simulated data. The results show that it allows a robust reconstruction of promoter activities, growth rates and protein concentrations from fluorescence and absorbance signals [Zulkower et al., 2015].

## 4.3 Discussion and perspectives

In this chapter, we have seen through two case studies how regularized linear inversion methods allow the reconstruction of interesting quantities from metabolomics or gene expression data. The two methods are able to capture the abrupt changes of metabolism and gene expression
during growth transition, and to avoid over-fitting the data when the system dynamics is slow. They are also able to deal with low signal-to-noise ratios such as the weakness of absorbance signal at the beginning of the experiment. Two different regularization methods haven been chosen: a stochastic approach based on a Bayesian framework for the metabolomics data and a deterministic approach using Tikhonov regularization for the reporter gene data. The approaches can be extended to retrieve time-varying profiles in a much wider range of problems, such as single-cell measurements of fluorescence signals. A necessary condition for their application is that the data should depend linearly on the signal to be estimated.

An important aspect of these approaches is that they come with a software implementation [Cinquemani et al., 2017, Martin et al., 2019, Zulkower et al., 2015], which makes them accessible to a larger audience. Eugenio Cinquemani and I are considering the development of new versions of the software for the analysis of metabolomics data. We envision to enhance the robustness of the method to noise in data at the beginning of the growth kinetics, and to include a graphical user interface for use by non-experts. This will be the object of an Inria technological development project proposal to be submitted in the near future.

Data analysed with these approaches were crucial for the work described in the rest of the manuscript. The following chapter provides one such example. Reporter gene data were used to understand cell physiology and based on the insights this provided, we proposed a strategy to control the bacterial growth.
<figure><img src="image_18.png" /></figure>
# Analysing and controlling cell physiology

The work presented in this chapter is the unexpected development of the analysis of the carbon starvation response described in Chapter 3. Follow-up studies made us realise the major role played by the global control of gene expression in the adaptation of bacterial growth, resulting from adjustments in the activity of the gene expression machinery. This work on global control was performed as part of the PhD thesis of Sara Berthoumieux. I co-supervised the early work of Sara on the global control of gene expression during her Master thesis and was closely involved in this part of her PhD project for the analysis and interpretation of the data. I also designed and performed some of the genetic constructions used in the study. This work and other work published at the same time, drew attention to the fact that the functioning of biochemical networks cannot be disconnected from the physiological state of the cell [Berthoumieux et al., 2013].

The observation motivated the follow-up study and the bet that adjusting the intracellular level of either ribosomes or RNA polymerase could be used to control the growth rate. This work was performed by Jérôme Izard as part of his PhD thesis, whom I co-supervised with Hans Geiselmann and Stéphan Lacour, and by the post-doctoral researcher Cindy Gomez Balderaz, co-supervised by Hans Geiselmann and Hidde de Jong. While controlling the synthesis of ribosomal proteins was not met with success, the strategy to control expression of the two limiting $\beta\beta'$ subunits proved more promising [Izard et al., 2015]. In addition to the co-advisorship and research design, I performed some experiments myself and contributed to the data analysis, notably the microfluidics data. This study was my first encounter with synthetic biology approaches. It was also a long and rich scientific adventure, which continues nowadays through a new line of research on our team agenda, dedicated to the analysis of resource allocation problems (see Chapter 8).

## 5.1 Contribution of cell physiology to the global control of gene expression

In the work discussed in Chapter 3, we modelled the various genetic regulatory interactions involved in the response of *E. coli* to carbon starvation. The model summarized the view that gene expression changes accompanying adaptation of bacterial growth to the stress result from combined effects of positive and negative transcription regulators. In the absence of quantitative and dynamical data in the literature that could be confronted to the model, J. Geiselmann and his group launched in the wake of this study a new experimental program, with my help. It aimed at quantifying the dynamics of gene expression in *E. coli* by means of reporter genes. Using the approach described in Chapter 4, we reconstructed promoter
activities from the fluorescence and absorbance measurements. An example is provided in Figure 5.1 for genes of a subpart of the carbon starvation response network. They include the global regulators of transcription Fis and CRP, as well as the metabolic gene *acs* coding for the acetylCoA synthetase. This enzyme allows the use of acetate as a carbon source. It is inhibited by Fis and activated by the complex CRPcAMP. Our former study suggested that Fis and CRP are part of a regulatory switch controlling growth adaptation, with each gene inhibiting each other's expression, allowing the expression of Fis in exponential phase and that of CRP in stationary phase. The promoter activities in Figure 5.1 tell a different story. The profiles for CRP and Fis are similar, with high levels in exponential phase that drop as cell growth rate slows down. The observation of an induction of *Acs* at the entrance in stationary phase is consistent with other reports [Wolfe, 2005]. We also obtained a similar gene expression profile for Fis in a former study with a different experimental approach [de Jong et al., 2010].

The similarity of gene expression profiles for *fis* and *crp* suggests that a common control mechanism might be at work. The two genes share at least the same gene expression machinery, whose activity is known to vary with the physiological state of the cell, notably the growth rate in steady-state conditions [Dennis and Bremer, 2008]. The concentration of available RNA polymerase is one example, but this quantity is difficult to monitor directly [Klumpp and Hwa, 2008]. As an indirect read-out of the global physiological state we therefore decided to use a constitutive promoter, the pRM promoter of phage λ, whose activity is controlled by the transcription and translation machinery and the pools of precursor metabolites, but not by any particular transcription factor in *E. coli* [Berthoumieux et al., 2013]. The activity of the promoter is shown in Figure 5.1(e). It is stationary in exponential phase and stabilizes to a lower level when cells cease growing. Considering this profile as representative of the physiological state of the cell, we developed an approach to dissect the contribution of the transcription factors and the global physiological state to gene expression. The promoter activity was then formulated as follows:

$$p(t) = k p_1(t) p_2(t), \qquad (5.1)$$

with $k$ [M min⁻¹] representing the maximum promoter activity. The dimensionless term $p_1(t)$, for convenience assumed to vary between 0 and 1, quantifies the modulation of the promoter activity by the global physiological state, for instance through the availability of free RNA polymerase. The dimensionless term $p_2(t)$, also varying between 0 and 1, accounts for the effect of transcription factors and other specific regulators, and may take the form of sigmoidal regulation functions as seen in Chapter 3. Normalization of the promoter activities with a reference state $t^0$ - chosen here at growth transition - and a subsequent logarithmic transformation gave the following expression:

$$\log \frac{p(t)}{p^0} = \log \frac{p_1(t)}{p_1^0} + \log \frac{p_2(t)}{p_2^0}, \qquad (5.2)$$

with $p_1^0 = p_1(t^0)$ and $p_2^0 = p_2(t^0)$.

This simple model can be used to test different hypotheses on the contribution of global and specific effects on gene expression. For instance, dominance of the global physiological effect makes the effect of transcription factors negligible. This is translated by $p_2(t) \approx p_2^0$ and thus
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

**Figure 5.1** – (a) Central regulatory circuit involved in the control of *E. coli* carbon metabolism, consisting of the two pleiotropic transcription factors Crp and Fis and their regulatory interactions. The global physiological state affects the expression of all genes in the network. Genes are shown in blue and promoters in red. Specific regulatory interactions are indicated by dashed lines and the effect of the physiological state by solid lines. (b-e) Experimental monitoring of transcriptional response of network. Time-varying activity of *fis* promoter (b, blue), derived from GFP reporter data, and absorbance (solid line, red). (c-e) Idem for the activities of the *crp* and *acs* promoters, as well as the activity of the pRM promoter of phage λ. The latter promoter is constitutive in the conditions studied and reflects the global physiological state of the cell. From [Berthoumieux et al., 2013].
$\log(p_2(t)/p_2^0) \approx 0$ in the right-hand side of Equation (5.2). In Figure 5.2(a-c), for instance, we test the contribution of global physiological effects (represented by pRM activity) to the promoter activity of *fis*, *crp*, and *acs* from Figure 5.1. Variation of the log normalized activity of pRM (representing the contribution of global physiological effects) explains the variation of the log normalized activity of promoter *fis* and *crp*, while it does not for *acs*. Taking into account both global and specific effects are necessary for this gene. It is induced at the onset of stationary phase, when its transcription factor CRP is itself activated by cAMP. Using dynamic data on intracellular cAMP concentration obtained from measurements of external cAMP (Figure 5.2(e)), we were able to show that the log normalized variation of intracellular cAMP concentration explains the variation of the remaining log normalized activity of *acs* promoter, once we have subtracted the contribution of the global physiological effect. These results and others obtained with variant strains in [Berthoumieux et al., 2013] show the importance of taking into account global physiological effects in gene expression. The specific regulatory effects by transcription factors appear to fine-tune the gene expression level set by global effects. This new vision on gene expression was confirmed in parallel for different promoters in different environmental conditions and micro-organisms (*E. coli* and yeast) [Gerosa et al., 2013, Keren et al., 2013]. These various works have revived past work in the 1950s-70s on the dependency of the macromolecular composition of the cell on growth rate, in particular the RNA and protein content [for review Dennis and Bremer, 2008]. They also paved the way to a synthetic biology approach in a follow-up study, aiming at controlling bacterial growth.

## 5.2 A synthetic biology approach to control bacterial growth

In the previous section, we analysed the natural strategies evolved by bacteria to adjust gene expression to growth changes through a modification of the activity of the gene expression machinery. This allows reallocating the available resources to the production of RNAs and proteins in the proportions needed at the new growth rate. Changing these resource allocation strategies has many applications in both fundamental microbiology and biotechnology. As we reviewed in [de Jong et al., 2017b], some of them re-engineer the transcription and translation machinery to re-allocate resources towards the production of proteins or compounds of interest. In [Izard et al., 2015], we proposed such an approach, by modifying the natural control of the synthesis of the two limiting RNA polymerase $\beta\beta'$ subunits in *E. coli*.

We engineered an *E. coli* strain, in which we replaced the natural promoter of genes *rpoBC* by a synthetic *lac* promoter inducible by IPTG and decoupled the synthesis of ribosomal proteins encoded by genes *rplKAJL* from that of the $\beta\beta'$ subunits (Figure 5.3(a); see [Izard et al., 2015] for more details). We also added two chromosomal copies of the *lacI* gene to confer mutational robustness to the strain. Inhibition of the synthetic *lac* promoter exerted by the *lac* repressor is relieved by IPTG addition. This leads to the expression of the two limiting $\beta\beta'$ subunits allowing RNA polymerase formation and cell growth. We characterized the physiology of the strain in [Izard et al., 2015] and the main results are shown in Figure 5.3. At high IPTG concentrations, the growth curves of the engineered "R" strain and the reference wild-type "W" strain are comparable (Panel (b)), as well as their growth rates at steady state in exponential phase (Panel (c)). The reference strain "W" is a wild-type strain of *E. coli* including the two extra copies of *lacI* gene. At low IPTG levels, cells stop growing after some time and their
<figure><img src="image_24.png" /><figcaption>**Figure 5.2** – Predicted and observed control of *fis*, *crp*, and *acs* activity by CRPCAMP and the physiological state of the cell, in various experimental conditions and genetic backgrounds. (a) Predicted (–, black) and measured (•, blue) relative activity of the *fis* promoter (log(*p*<sub>fis</sub>(*t*)/p<sup>0</sup><sub>fis</sub>)) as a function of the relative activity of the pRM promoter (log(*p*<sub>RM</sub>(*t*)/p<sup>0</sup><sub>RM</sub>)). (b-c) Idem for *crp* and *acs*. (d) Predicted (–, black) and measured (•, blue) remaining relative activity of the *acs* promoter after subtraction of the effect of global physiological parameters (log(*p*<sub>acs</sub>(*t*)/p<sup>0</sup><sub>acs</sub>) – log(*p*<sub>RM</sub>(*t*)/p<sup>0</sup><sub>RM</sub>)) and as a function of the relative intracellular cAMP concentration (log(*c*(*t*)/c<sup>0</sup>)). (e) Absorbance (red) and derived concentration of intracellular cAMP from measurements of the external cAMP concentration. The confidence intervals in the plots have been computed from experimental replicas. From [Berthoumieux et al., 2013].</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 5.3** – (a) Construction of an *E. coli* strain with inducible expression of the *rpoBC* genes encoding the *ββ'* subunits of RNA polymerase. (b) Growth kinetics in a microplate of the engineered "R" strain (blue) and the wild-type strain "W" (red) at different concentrations of IPTG added to M9 minimal medium supplemented with 0.2% glucose. The growth rates are typically computed in the time intervals indicated by the green bars. (c) Growth rate estimated from the absorbance data for the R and W strains. The black circles represent the growth rates obtained with the R and W strains in shake flask experiments. (d) Quantification of the *ββ'* subunits of RNA polymerase using a chromosomal fusion of the *rpoC* gene with the gene encoding the mCherry fluorescent reporter protein. (e) Quantitative dependence of the growth rate on *ββ'* concentrations, computed from the data in (d) and the measured growth rates at the different IPTG concentrations. The blue curve is a Hill function with a Hill coefficient of 10. From [Izard et al., 2015].</figcaption></figure>
growth rate is null. These results obtained in M9 minimal medium supplemented with glucose were also observed in other growth media and genetic backgrounds. Labelling protein $\beta'$ in the R and W strains with the fluorescent protein mCherry allows quantifying the expression of the RNA polymerase subunits (Panel (d)). A high IPTG concentration results in the overexpression of the subunits compared to the wild-type situation. At low IPTG, the subunit levels is too low to support growth (Panel (e)). We have thus a growth switch, with growth that can be turned on or off by the simple addition or removal of inducer.

Time-lapse microscopy in a microfluidics experiment followed by the quantification of individual cell growth rates showed that the switch is reversible and that almost all growth-arrested cells resume growth when IPTG is added back to their medium (Figure 5.4(a)). Using the engineered strain allows to reallocate resources from biomass formation to the production of a metabolite of interest, glycerol in this case. We transformed the R and W strains with a plasmid carrying the genes that code for the glycerol pathway in yeast. Following the growth arrest and the depletion of RNA polymerase, the production yield of glycerol in the R strain is twice higher than in the W strain (Figure 5.4(b)-(d)), and close to the maximal theoretical yield [Izard et al., 2015]. This example shows that growth-arrested cells remain metabolically active and represent an interesting platform for biotechnological applications. In a follow-up study, we showed in collaboration with the company Metabolic Explorer that the synthetic growth switch still works efficiently in liter-scale bioreactors (Ropers et al., submitted). We also characterized the physiology of the growth-controlled strain, by assessing the reorganization of the metabolome, transcriptome, and proteome induced by the growth switch. The data confirms the idea that the growth-controlled strain is functioning like a bag of active metabolic enzymes.

## 5.3 Discussion and perspectives

The global control mechanism evidenced in our study in [Berthoumieux et al., 2013] is the consequence of changes in the pools of ribosomes and RNA polymerases, which affects transcription and translation rates directly, modifies gene expression on a global scale, and ultimately affects cell macromolecular composition (reviewed in [Jun et al., 2018]). Our work shows that these observations, originally made for *E. coli* cells in steady-state growth, remain valid in dynamic environments: when growth slows down following glucose exhaustion in the growth medium, changes in the expression of global regulators of transcription can be attributed to these global effects rather than to specific regulations [Berthoumieux et al., 2013]. This shows that the functioning of biochemical networks cannot be disconnected from the physiological state of the cell [Berthoumieux et al., 2013]. Following work took this phenomenon into consideration to analyse the cellular metabolism and mRNA decay (see Chapters 6 and 7).

The exploitation of global control of gene expression for biotechnological purposes allowed us to engineer an *E. coli* strain and turn it into a growth switch [Izard et al., 2015]. I am still full of wonder at the results obtained with this strain. While flexibility characterizes the metabolism of bacterial strains and often oppose modifications aiming at hijacking it for bio-production purposes for instance, the gene expression machinery resembles an Achille's heel. It has such a profound impact on cell processes, that it becomes possible to intervene directly in the resource
<figure><img src="image_26.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_30.png" /></figure>

**Figure 5.4** – (a) Growth arrest by external control of *rpoBC* is reversible. The R strain was grown in a microfluidics device and phase-contrast images were acquired every 10 min. Cells were grown in M9 minimal medium with glucose, initially in the presence of 1M IPTG. 6 h after removing IPTG from the medium, the cells are elongated. About 100 min after adding back 1M IPTG into the medium, the elongated cells divide and resume normal growth. The growth rates in the plot are the weighted mean of the growth rates of 100 individual cells. The glycerol-producing W strain (W-gly) (b) and R-gly strain (c) were grown in shake flasks in M9 minimal medium supplemented with 2 g/L glucose. The optical density (OD<sub>600</sub>, black), the glycerol (red), and glucose (blue) concentrations were measured in samples taken at intervals of about 30 min using coupled enzyme assays. (d) The instantaneous yield of glycerol production in the W-gly is computed by dividing the glucose consumption rate by the glycerol production rate, in a time interval in which the derivatives of the glucose and glycerol concentration curves are well defined (between 250 and 450 min). (e) Idem for the R-gly strain (in the time interval between 250 and 700 min).
allocation scheme by targeting the machinery as we did in our study. The engineered strain is a great tool, which has utility both in fundamental research and biotechnology. On the one hand, it allows the analysis of the natural resource allocation strategies that have evolved and, on the other hand, it constitutes a new method for producing metabolites, peptides and recombinant proteins, which we patented in Europe and the USA [Geiselmann et al., 2015].

The work described so far focused on small subnetworks controlling the adaptation of *E. coli* to a carbon source. The main reason is that we restricted our analysis to measurable parts of the *E. coli* network, for which data was hence available. However, the increased availability of omics data makes it possible to analyse biological processes involved in bacterial adaptation at a larger scale. This will be the subject of the two following chapters dedicated to the analysis of the regulation of metabolism and mRNA decay at the genome-wide level in *E. coli*.
<figure><img src="image_31.png" /></figure>
# Metabolic network models as platforms for integrating omics data

This chapter is dedicated to the analysis of the regulation of metabolism in *E. coli* using constraint-based models, which augment the information content of the data. The work presented in this chapter was performed as part of the PhD thesis of Manon Morin funded by an INRA-Inria PhD grant. Manon was co-supervised by Muriel Cocaign-Bousquet, Brice Enjalbert, and myself. Muriel was the main supervisor and I lead and performed the modelling work and data integration shown in the two studies. With this work, I made my first steps in constraint-based modelling and the exploitation of metabolomics data. This is the reason why I chose to highlight the two corresponding papers [Morin et al., 2016, 2017]. I remain impressed to see how much information can be extracted from the data, based on the sole information on the measured uptake and secretion rates, as well as the growth rate, plus additional intracellular metabolite concentrations. I used a similar approach in the framework of the PhD thesis of Stéphane Pinhal, whom I co-supervised with Hidde de Jong and Johannes Geiselmann, dedicated to the analysis of growth inhibition by acetate [Pinhal et al., 2019].

## 6.1 Post-transcriptional regulation of central carbon metabolism in *E. coli*

The possible influence of mRNA stability on metabolic activity has long been ignored. With the identification of global regulators such as the carbon storage regulator system (CSR) and of small RNAs controlling the stability and/or translation of metabolic genes, it has become clear that post-transcriptional regulations add a new layer to the already complex network controlling cell metabolism in *E. coli* [Kotte et al., 2010, Kochanowski et al., 2013]. To what extent this added complexity is determinant in the adaptation of bacterial growth to the environment was unclear when we started this study. We focused on the CSR system, whose functioning is paradoxical. CSR consists of the dimeric mRNA binding protein CsrA and small regulatory RNAs CsrB/C, which inhibit CsrA activity. Both of these noncoding RNAs are targeted by the protein CsrD, which triggers their RNase E-dependent degradation (see Pourciau et al. 2020 for a recent review). The system is essential for growth on glycolytic media, where pleiotropic regulators like ppGpp, CRP, and RpoS have larger regulons but are not essential. Past studies conducted with a multitude of different strains and conditions showed that attenuating the *csrA* gene through deletion of the last 10 amino acids leaves cells viable, albeit with a strongly reduced growth and perturbed biofilm formation, motility, and the accumulation of the storage sugar glycogen [Esquerré et al., 2016, Romeo et al., 1993]. Our objective was to obtain a more global view of the CSR regulon and understand the physiology of the attenuated *csrA* strain and other CSR deletion mutants. Is really the accumulation of glycogen responsible for the
essentiaity of CsrA, as proposed in the literature [Timmermans and Van Melderen, 2009]? Or
is there an alternative explanation to this phenomenon? Do post-transcriptional regulation play
an important role in shaping the central carbon metabolism? We characterized CSR targets
among the central carbon metabolism genes to obtain responses to these questions, through a
multi-scale analysis of growth properties, mRNA levels, enzyme activities, fluxes and metabolite
concentrations in CSR variant and wild-type strains, as summarized in Figure 6.1.

<figure><img src="image_32.png" /><figcaption>**Figure 6.1** – Systematic measurement of different molecular levels in wild-type and mutant/deletion strains of the CSR system during a glucose-acetate diauxie in *E. coli*.</figcaption></figure>

In a first study, we probed the activity of central carbon metabolism in the attenuated *csrA* and wild-type strains growing exponentially in minimal M9 medium supplemented with glucose [Morin et al., 2016]. Figure 6.2 shows in bold the names of metabolites and reactions in glycolysis/gluconeogenesis, glycogenesis, pentose-phosphate pathway, and the Krebs cycle for which experimental data was obtained in the study - a measured concentration, enzyme specific activity and/or mRNA level obtained by qRT-PCR. The corresponding data set is shown in Figure 6.3. The flux through each reaction was determined by metabolic flux analysis. In the absence of transcriptomics data at the time (they were obtained later on, [Morin et al., 2020]), we could not develop a condition-specific model from a generic metabolic reconstruction [Bordbar et al., 2014]. Instead we tried to exploit as much as we could the available metabolomics data to develop models as specific to the strains and conditions studied as possible. We adapted the iAF1260 genome-scale reconstruction of *E. coli* metabolism [Orth et al., 2010] for the metabolism and storage of glycogen, without distinguishing between the different forms of the molecule, because only the total pool of glycogen was quantified. Maintenance fluxes were determined from previous measurements of steady-state growth rates and glucose uptake rates [Esquerré et al., 2014].

We then formulated a limited number of uptake and secretion constraints that were directly
motivated by the composition of the growth medium and the utilization of glucose as the sole
carbon source and the secretion of acetate and CO₂. These were determined by means of the approach described in Chapter 4. Intracellular metabolite concentrations were used to specify thermodynamic constraints that enforce reaction directionality, based on the determination of the Gibbs energy of the reactions at room temperature, and physiological pH and ionic strength [Fleming et al., 2009, Henry et al., 2006, 2007]. The resulting constraint model was then used to predict the intracellular metabolic fluxes. Instead of using a principle of optimality to solve the model, we performed a metabolic flux analysis, where we sought to minimize the difference between the measured and predicted fluxes, including the biomass production rate. This was formulated as a linear programming problem, similarly to what was done in [Lee et al., 2012]:

$$
\begin{align*}
\min \quad & \sum_{j=1}^{p} (u_j^+ + u_j^-) \text{ subject to:} \\
& N v = 0, \\
& v^l \le v \le v^u, \\
& v_j - u_j^+ + u_j^- = u_j^M, \quad \text{for all } j = 1, \dots, p, u_j^+, u_j^- \ge 0. \quad (6.1)
\end{align*}
$$

with $v$ the vector of steady-state fluxes with lower bounds $v_l$ and upper bounds $v_u$, $N$ the stoichiometry matrix, $u_j^+$ and $u_j^-$ non-negative dummy fluxes, and $u_M$ the vector of $p$ measurements of exchange fluxes and growth rate. We assumed that the first $p$ elements of $v$ correspond to the measured fluxes. We further analysed the optimal solutions by flux variability analysis [Mahadevan and Schilling, 2003] to determine the minimum and maximum flux values satisfying the constraints and consistent with the measurements. This allowed restricting the possible values of intracellular fluxes to very tight intervals, as shown in Figure 6.3.

A genetic perturbation such as CsrA attenuation may affect metabolism in different manners: directly, through the modification of the expression of CsrA target genes, or indirectly through modifications of the growth rate or of upstream or downstream reactions, which are propagated through the metabolic pathway (Figure 6.3(e)). The effect of growth rate on metabolism could be assessed by using data from a previous study on continuous cultures growing at rates comparable to the attenuated and wild-type strains (0.3 and 0.6 h⁻¹) [Esquerré et al., 2014] (see Figure 6.3). We separated the direct and indirect effect of CsrA attenuation on metabolic fluxes using a hierarchical regulation analysis [ter Kuile and Westerhoff, 2001, van Eunen et al., 2011]. The approach quantifies the contribution of changes in gene expression or metabolite pool to the flux change. We write the rate of an enzyme-catalysed reaction at steady-state $J = v(e, x, K) = f(e) \times g(x, K)$, in which $v$ is the rate, $e$ the enzyme concentration, $x$ the vector of metabolite concentrations (substrates, products, and effectors), and $K$, the vector of dissociation constants. A logarithmic transformation dissects the flux into two terms, one depending on the enzyme concentration and the other, on the metabolite concentrations: $log J = log f(e) + log g(x, K)$. The change of flux between the wild-type and *csrA51* strains is written: $\Delta log J = \Delta log f(e) + \Delta log g(x, K)$. The contribution of gene expression and metabolic control to the flux change is quantified by dividing the latter expression as follows:

$$ \frac{\Delta \log J}{\Delta \log J} = \frac{\Delta \log f(e)}{\Delta \log J} + \frac{\Delta \log g(x, K)}{\Delta \log J} = \rho_h + \rho_m = 1, \quad (6.2) $$
in which $\rho_h$ is the hierarchical regulation coefficient and $\rho_m$, the metabolic regulation co-efficient. In our conditions, $f(e)$ corresponds to the measured specific activity SA of the enzyme. We can therefore determine $\rho_h$ directly from the experimental data and the fluxes: $\rho_h = \Delta\log SA/\Delta\log J$, and deduce $\rho_m = 1 - \rho_h$. For each reaction, sets of hierarchical regulation coefficients were calculated from the specific activities of the different replicates of each strain, and the lower and upper bounds for the fluxes, which gives the boxplots in Figure 6.3(f).

These experimental and model analyses allowed us to demonstrate the strong control of the upper part of glycolysis by the CSR post-transcriptional regulatory system [Morin et al., 2016]. Attenuation of CsrA activity results in a decrease of most glycolytic activities, especially the phosphofructokinase. This causes an accumulation of metabolites in the upper part of glycolysis before the phosphofructokinase step and results in a glucose-phosphate stress controlling negatively the sugar uptake. This strongly affects the bacterial growth rate and could explain the essentiality of *csrA* gene for growth on glycolytic substrates. The glucose-phosphate stress can be relieved by restoring PfkA activity in the *csrA* mutant strain.

## 6.2 Post-transcriptional regulation of metabolic adaptation

Our previous study demonstrated the major role played by the CSR system in the control of upper glycolysis, but its putative role during metabolic adaptation remained to be established. In addition we showed that metabolic reprogramming rather than glycogen accumulation causes a growth defect when CsrA is attenuated: does glycogen have a physiological role in *E. coli*? This polysaccharide is the main storage form of glucose, from bacteria such as *Escherichia coli* to yeasts and mammals. Although its function as a sugar reserve in mammals is well documented, the role of glycogen in bacteria is not that clear. In a follow-up study using the same data set obtained with our multi-scale analysis of central carbon metabolism (Figure 6.1), we analysed gene expression and metabolic pools in CSR variants and wild-type strain transitioning from growth on glucose to growth on acetate [Morin et al., 2017].

The main results of this study are summarized in Figure 6.4. The figure illustrates the evolution of biomass and pools of glucose, acetate and glycogen, as cells grow first on glucose and then transition to acetate when the preferred carbon source is depleted (Panels (a)-(d)). Glycogen is also consumed during the second growth phase. The timing of acetate consumption varies between strains (Panel (c)). We could not attribute it to differences in gene expression levels, but to differences in the energetic status. Analyses of the adenylate energy charge $AEC = [ATP + 1/2 ADP]/[ATP + ADP + AMP]$ before glucose depletion showed a high energetic status, typical of exponential growth, while it drops after glucose depletion, to low levels in the *csrBC* strain reminiscent of dying cells (Panel (e)). Using the metabolic model previously developed in [Morin et al., 2016], we predicted by flux balance analysis the maximal ATP fluxes that the various wild-type and CSR variant strains are able to produce following glucose depletion, in the range of measured glycogen and acetate rates (Panel (f)). Acetate and glycogen appear to cover the needs for maintenance energy. Preventing glycogen synthesis through deletion of the *glgC* gene does not allow cells to maintain their energetic status after glucose depletion (Panel (g)). This study brings a new vision to the physiological role of glycogen. The stored polysaccharide
<figure><img src="image_33.png" /><figcaption>**Figure 6.2** – Schematic representation of the main reactions in central carbon metabolism including the glycolysis, pentose-phosphate pathway, glycogenesis, and tricarboxylic acid (TCA) cycle. Metabolites in bold were quantified by metabolomics. Reactions shown in bold were monitored, by qRT-PCR measurement of the mRNA levels and/or determination of the enzymatic activities. Metabolic network analysis allowed to predict the flux value through each metabolic reaction. Dashed lines represent aggregated reactions. 2/3PG corresponds to the pool of 2PG and 3PG that cannot be distinguished in metabolomics experiments. The corresponding data is shown in Figure 6.2.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Figure 6.3** – Comparison of the (a) fluxes, (b) mRNA levels, (c) enzyme activities, and (d) metabolite concentrations in the central carbon metabolism between the wild-type and the *csrA51* strain or at different growth rates. The displayed values correspond to the log2 of the ratio of the *csrA51* strain to its isogenic wild type. Hatched columns represent log2 of ratio between values in the wild-type growing at $\mu = 0.3h^{-1}$ (chemostatic cultures) to the same strain at $\mu = 0.6h^{-1}$ (batch cultures) (data from Esquerré et al. 2014). (e) Schematic representation of direct and indirect effects of CsrA on metabolic fluxes. (f) Hierarchical and metabolic regulation coefficients are shown as boxplots. Boxes represent the interquartile range (IQR) between the first and third quartiles. Whiskers denote the lowest and highest values within 1.5 IQR from the first and third quartiles. Purple, median metabolic coefficient; orange, median hierarchical coefficient.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Figure 6.4** – Behaviour of CSR system mutants during the glucose-acetate transition. Changes in (a) the biomass concentration, (b) extracellular glucose concentration, (c) extracellular acetate concentration and (d) glycogen concentration. All data concerning the replicates are displayed as dots, and the fitted average value of each strain is displayed as a line. Shaded areas represent ±1 standard deviation. WT, black circles; *csrBC* mutant, green squares; *csrD* mutant, orange triangles; *csrA51* mutant, red diamonds. (e) Adenylate energy charges (AEC) of the four strains, during glucose consumption (plain bars) and 1.5 h after glucose exhaustion (striped bars). A significant difference between a mutant and the WT is represented by an asterisk (P < 0.05 [t-test]). (f) Maximum flux of ATP predicted by the constraint-based model in [Morin et al., 2016]. The bar to the right of the plot indicates the ATP flux values. The four strains were assessed at two different times 90 min (*) and 150 min (**). (g) Effect of the deletion of *glgC* gene on the AEC of the wild-type and *csrA51* strains.</figcaption></figure>
is utilized during transition to new carbon source to provide energy for growth resumption. As
major regulator of glycogen storage, the CSR system appears to play a key role for the fitness
of *E. coli* cells transitioning from growth on glycolytic to gluconeogenic substrates.

**6.3 Discussion and perspectives**

Integration of data with models augments the information content of the data. The two
studies presented here provide one such example: the constraint-based model integrating
metabolomics measurements allowed us to determine intracellular fluxes, while small models in
hierarchical regulation analyses integrating the fluxes and enzyme specific activities enabled us
to disentangle the direct and indirect regulatory effects of CsrA.

It is, however, difficult to causally relate the predicted changes of fluxes to changes in
gene expression, metabolism, and other interesting parameters such as growth rate. While the
approach in [Morin et al., 2016] enabled us to perform this analysis for a subset of reactions,
the next step will be to develop an approach allowing to identify at the genome-scale level the
specific parts of metabolism directly affected by a perturbation. I will come back to this point
in Chapter 8, where I will discuss some work started on this subject.

The field of constraint-based modelling is evolving fast, notably with the development of
approaches allowing the integration of omics data [Bordbar et al., 2014]. The recent obtention
of transcriptomics data such as those obtained in the same experimental conditions [Morin et al.,
2020] will help us improve our analysis of the post-transcriptional regulations of cell metabolism.
Work in this direction has recently started with the internship of Amélie Caddéo (Univ. Grenoble
Alpes), whom I am supervising. I will discuss this specific point in Chapter 8, for a new project
in which we seek to understand the cross-talk between RNA metabolism and central carbon
metabolism.
# Analysis of bacterial mRNA decay

This chapter is devoted to the study of mRNA degradation kinetics. Transcriptomics experiments by means of microarrays or RNAseq enable the monitoring of the decay of cellular mRNAs at the genome scale following a transcription arrest. Simple models assuming exponential decay are generally used to determine mRNA half-lives from the degradation profiles [Laguerre et al., 2018]. However, the assumption that mRNA decay follows first-order kinetics does not permit a detailed investigation of the regulatory mechanisms responsible for the diversity of degradation profiles observed experimentally. Can we develop a model of mRNA degradation that would be a good compromise between a mechanistic description of the degradation process and simplicity for the confrontation of the model with experimental data?

We addressed the question in the framework of the PhD thesis of Thibault Etienne, funded by an INRA-Inria PhD grant, which I co-advised with Muriel Cocaign-Bousquet. I was the main supervisor of Thibault. We made the bet that more information on the regulatory mechanisms could be extracted from dynamic omics data by confronting them to mechanistic models of mRNA degradation instead of exponential models. Our final objective was to be able to obtain more information than just the determination and analysis of mRNA half lives from a large data set obtained in a former study monitoring the degradation of 4254 cell mRNAs in steady-state cultures of *E. coli* growing at four different rates [Esquerré et al., 2014]. The corresponding work is described below. In a first section, I describe the mechanistic modelling of mRNA degradation and how we used the model to show that competition between mRNAs could affect the degradation kinetics [Etienne et al., 2020]. In a follow-up study described in Section 7.2, work that will be submitted soon, we show the physiological relevance of this phenomenon using dynamic omics data [Etienne et al., In preparation]. To the best of my knowledge, this work on mRNA decay is the first example of the interpretation of transcriptomics data by means of a mechanistic model that allows to identify regulatory mechanisms and analyse their contribution to the cell physiology.

## 7.1 Competitive effects in bacterial mRNA decay

One of the most common approaches for the experimental determination of mRNA half-life is the monitoring of residual mRNA concentrations following transcriptional arrest. An antibiotics like rifampicin is often used, which blocks the elongation of transcription by RNA polymerase. After a delay during which RNA polymerase completes the mRNAs it had started to transcribe before antibiotics addition, mRNAs no longer accumulate and are progressively degraded by the cell machinery called degradosome (Figure 7.1, yellow panel). Sequence and structure characteristics of the mRNAs are factors known to result in variations of the degradation rate, as well as regulations by small RNAs and/or RNA-binding proteins such as Hfq. Since many of the degradation profiles resemble an exponential decay, the curves are fitted by an exponential
function and the estimated rate constant is used for the calculation of mRNA half-lives. This method and some variants are the reference for the determination of mRNA half-life [Laguerre et al., 2018]. Alternative modelling approaches have analysed the coupling of degradation with other cellular processes, but they are too detailed to allow a thorough confrontation of model predictions to data (for review, Roux et al., submitted).

The development of models of mRNA degradation and their analysis is described in [Etienne et al., 2020]. Based on the experimental literature, we showed that we can assimilate mRNA decay by the degradosome to a macro-reaction catalysed by the major enzyme RNase E. While degradation of each mRNA is usually studied independently from other mRNAs, we realised that cellular mRNAs share the same machinery of degradation and should compete for binding to RNase E because they outnumber the enzymes within cells. Modelling of the degradation process with mass-action law led to two models, including the mechanism of competition or not. At the moment of reducing the models, we noticed that the standard QSS approximation was not applicable in our conditions. First, because the total and not the free mRNA concentration is most likely measured in the transcriptomics assays (mRNAs are deproteinized during cell extraction). Second, because available literature data indicate an excess of RNase E with respect to individual mRNAs, while the total concentration of all cell mRNAs is way larger than the enzyme concentration. These conditions motivated the use of the total QSSA described in Section 2.1.3, since it considers the total substrate concentration and has a larger domain of validity than sQSSA, which should facilitate the future step of parameter estimation. Application of the first-order tQSSA gives for the model describing single mRNA degradation, in the absence of competition:

$$ \frac{d}{dt} m_i(t) = - \frac{k_{cat} E_0 m_i(t)}{K m_i + E_0 + m_i(t)} \quad (7.1) $$

with, as seen before, the following domain of validity (Section 2.1.3):

$$ E_0 + K m_i \gg m_{i0} \quad \text{and} \quad K_i \ll K m_i \quad \text{or:} \quad (7.2) $$

$$ E_0 \gg m_{i0} \quad \text{and} \quad E_0 \gg K m_i \approx K_i. \quad (7.3) $$

Building upon previous works by Pedersen et al. [2007] and Tang and Riley [2013], we developed an approximate version of the tQSSA form which includes competition between all cellular mRNAs:

$$ \frac{d}{dt} m_i(t) = - \frac{k_{cat} E_0 m_i(t)}{K m_i \left( 1 + \sum_{j \neq i} \frac{m_j(t)}{K m_j} \right) + E_0 + m_i(t)} \quad (7.4) $$

with $j = 1 \cdots n, n = 4254$ mRNAs, and any of the following conditions implying the validity of the approximation:

$$ E_0 \ll K m_i^{app}(0) + m_{i0} \quad \text{and} \quad K \lesssim K m_i^{app}(0) + m_{i0}, $$

$$ K m_i \gg \sum_{i=1}^{n} m_{i0} \quad \text{and} \quad K \ll K m_i^{app}(0), $$

$$ K m_i \gg \sum_{i=1}^{n} m_{i0} \quad \text{and} \quad E_0 \gg K \approx K m_i^{app}(0), $$

$$ E_0 \gg K m_i^{app}(0) + m_{i0} \quad \text{and} \quad E_0 \gg K, $$
<figure><img src="image_36.png" /><figcaption>**Figure 7.1** – Analysis of the physiological control of bacterial mRNA decay. Yellow panel: the objective of the study is to formulate hypotheses on the regulatory cellular mechanisms that could explain the observed adjustments of mRNA degradation profiles in *E. coli* cells growing at four different rates, like competition between mRNAs, protection of mRNAs by elongating ribosomes, or specific regulations involving small RNAs and RNA-binding proteins (CSR, Hfq,...). Orange panel: Nonlinear-mixed effect (NLME) modelling of bacterial mRNA decay using the degradation model in (7.4). We consider that the parameter values of individual mRNAs (Km values and initial concentrations) are drawn from a distribution common to the population of mRNAs. The black line represents the degradation profile of an average mRNA characterized by mean Km and initial concentration values. Blue panel: the parameters of individual mRNAs allow to generate degradation profiles that fit well the data. Green panel: parameters with values different from the population mean indicate the possibility of underlying regulatory mechanisms. Further analyses (e.g. functional annotation, enrichment analysis) are carried out to propose hypotheses on the possible regulatory mechanisms at work. From Etienne et al., In preparation].</figcaption></figure>
where $K = \frac{k_{cat}}{\min\{k_{1,i}\}}$. The degradation introduces a coupling between mRNAs. For the experimental conditions described in Figure 7.1 (yellow panel) for instance, this brings the model with competition to a large system of 4254 ordinary differential equations for each of the four environmental conditions. The model differs from the one in (7.1) by the multiplication of enzyme affinity with a competition term depending on the concentration of the other cell mRNAs. This shows that competition between mRNAs decreases enzyme affinity globally for all mRNAs.

To test if competition between mRNAs has any impact in the kinetics of degradation, we numerically simulated the models by means of parameters and initial conditions from the literature satisfying the above conditions [Etienne et al., 2020]. The profiles generated with the competitive model vary from seemingly exponential to more linear profiles (Figure 7.2(a-d)). Such diversity is difficult to reproduce in the absence of competition and is observed experimentally [Etienne et al., 2020]. The explanation lies in the titration of RNase E at the beginning of the kinetics and as long as the total mRNA concentration remains sufficiently high, which delays the onset of degradation (Figure 7.2(e,f)). During this period of time, mRNA competition slows down the degradation rate. By means of rate response coefficients evaluating the sensitivity of the degradation rate to changes in mRNA concentrations, we further showed that 1) competition differentially affects the fate of mRNAs – it stabilizes mRNAs with low affinities and destabilizes those with high affinities – and 2) it explains the observed negative correlations between mRNA concentrations and half lives indicating that more abundant mRNAs are degraded more rapidly [Esquerré et al., 2014, Esquerré et al., 2015, Nouaille et al., 2017].

The lag effect caused by mRNA competition was a striking result. Residual transcription is generally associated with delayed degradation, in particular for long genes or genes towards the 3' end of operons [Chen et al., 2015]. There is no direct evidence of this phenomenon in the literature, through measurements of the free or bound concentration of RNase E, for instance. To test this possibility, we analysed dynamical transcriptomics data in *E. coli*. We estimated the delay before degradation for 3140 mRNAs of this data set [Esquerré et al., 2014], as well as the maximal time needed to transcribe each of them using an elongation rate constant experimentally determined in [Chen et al., 2015]. Among the 2454 mRNAs that are not immediately degraded after rifampicin addition, 51% of them have a delay before degradation larger than the time needed for transcription. The delay is even twice for 21% of them (Figure 7.3). While this quick analysis clearly underestimates the number of mRNAs for which transcription elongation takes a shorter time than the delay before degradation, it indicates that residual transcription is not the sole determinant for the delay in the data set studied. Competition between mRNAs could well be another one. In a follow-up study, we verified the reality of the phenomenon of competition, by using data obtained in [Esquerré et al., 2014].

## 7.2 Integrative analysis of mRNA degradation

Esquerré et al. [Esquerré et al., 2014] showed that mRNA stability is one mechanism used by *E. coli* bacteria to adjust gene expression to their growth rate. Based on the data and our models, can we make hypotheses on the regulatory mechanisms at work? Is competition one of them? The principle of this study, which we are finalizing, is sketched in Fig. 7.1. We have shown in our study that nonlinear mixed effects (NLME) modelling [Lavielle, 2014] can be used to infer the parameters of the degradation model in (7.4) from dynamical transcriptomics data obtained by microarrays in *E. coli* cells growing at four different growth rates [Esquerré et al., 2014]. This framework generally yields good estimation results [Gonzalez et al., 2013].

In the NLME framework, we consider that the parameter values of individual mRNAs are drawn from a distribution common to the population of mRNAs. Concretely, we describe the time-series data
<figure><img src="image_37.png" /><figcaption>**Figure 7.2** – Numerical simulation of mRNA degradation kinetics in isolated and competitive systems. Predicted profiles for the (a) competitive and (b) isolated systems. The profiles are normalized to their respective initial concentrations for the competitive (c) and isolated (d) systems. Predicted free RNase E concentrations for (e) the competitive system ($E_{free} = E_0 - \sum_{i=1}^{n} c_i(t)$) and (f) the isolated one ($E_{i,free} = E_0 - c_i(t)$). 4312 curves are displayed in this case, due to the lack of coupling between mRNAs. The colour bars on the right side represent the normalized gradient of Km values, on a scale from zero (the minimal Km value) to one (maximal value). From [Etienne et al., 2020].</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 7.3** – Role of residual transcription in the retardation of degradation. On log scales and for each mRNA with a non negligible delay, the duration of transcription is plotted versus the delay before degradation. mRNAs with a delay before degradation higher than the time needed for transcription are displayed in blue and red, otherwise. From [Etienne et al., 2020].</figcaption></figure>

with the following measurement model:

$$Y_{i,t_{j,\mu},\mu} = f(m_{i,\mu}(t_{j,\mu}), \Phi_{i,\mu}) + g(f(m_{i,\mu}(t_{i,j,\mu}), \Phi_{i,\mu}), \theta) \times \varepsilon_{i,j,\mu}. \quad (7.5)$$

Here, $Y_{i,t_{j,\mu},\mu}$ is the observed concentration of mRNA i at time $t_{j,\mu}$ in a given growth condition $\mu$, $\varepsilon_{i,j,\mu} \sim \mathcal{N}(0, 1)$ is the residual error, and the function $g$ the residual error model with a vector of noise parameters $\theta$. Function $f$ is the solution to the ODE system (7.4). The vector of parameters $\Phi_{i,\mu}$ is a function of a vector of population parameters $\Phi_{pop}$, fixed effects $\beta_{\mu}$, and two vectors of random effects describing, respectively, the individual variability of mRNAs, $\eta_i \sim \mathcal{N}(0, \Sigma)$, resulting from differences in mRNA characteristics and regulations, and additional individual variability between growth conditions, $\eta_{i,\mu} \sim \mathcal{N}(0, \Omega)$:

$$\Phi_{i,\mu} = \Phi_{pop} + \beta_{\mu} \times \Phi_{pop} + \eta_{i} + \eta_{i,\mu}. \quad (7.6)$$

The estimation problem amounts to inferring the fixed effects ($\beta_{\mu}$) and parameter distributions describing the population (via parameters $\Phi_{pop}$, $\Sigma$, and $\Omega$), from which, together with the data, are then derived the specific parameters for individual mRNAs. This approach allowed us to estimate the individual kinetic parameters of the 4254 E. coli mRNAs in the four growth conditions. The physiological interpretation of the parameters is under way. Current results indicate that competition between mRNAs for binding to RNase E is indeed a global regulatory mechanism adjusting mRNA stability to growth rate. Additional specific regulatory mechanisms mediated by small RNAs or RNA-binding proteins (HFQ, CSR system...) appear to fine tune the stability of no more than a fifth of the mRNAs [Etienne et al., In preparation].

## 7.3 Discussion and perspectives

The modelling of mRNA decay described in this chapter has allowed the identification of a new regulatory mechanism of mRNA degradation. It relies on the competition between mRNAs for their binding to RNase E. In the first minutes after the onset of degradation, the enzyme is titrated by the myriad of cellular mRNAs. Sole those with higher concentrations and/or higher affinities have more chances to be degraded. Competition on the one hand tends to stabilize mRNAs by increasing the
competitive effect at the whole-cell level. On the other hand, it affects mRNAs individually through the modification of their affinity: competition stabilizes mRNAs with low affinities and destabilizes those with high affinities. This global mechanism allows to adjust degradation rates to the intracellular mRNA concentrations. The latter are in large part regulated at the transcriptional level [Esquerré et al., 2014, Morin et al., 2020]. This means that degradation is directly coupled to transcription and that mRNA competition is the cornerstone of this regulatory mechanism. Surprisingly, less than 20% of cellular mRNAs are the target of additional specific regulations.

The observation that the competition between mRNAs is a major global control mechanism, while specific regulatory mechanisms concern a small fraction of cellular mRNAs, is an unexpected result. This is nevertheless reminiscent of the global changes of the gene expression machinery studied in Chapter 5 that contribute more to transcription than specific regulatory mechanisms. It provides another example that biological processes cannot be studied without considering the physiological state of the cell. This work on mRNA decay was recently highlighted in a press release for the general public by INRAE<fnref n="1" />.

The NLME framework was critical in this study to allow the estimation of model parameters from microarray data and to cope with the noise in the data. While sub-sets of time-series data are generally analysed in studies of mRNA degradation, our approach is exhaustive since all cell mRNAs are included and their degradation profiles fitted. This allows to draw conclusions valid at the entire cell level. In addition, our approach is modular and scalable. It can be extended to other types of models and to other single-cell or population data, as well as additional heterogeneous data. For instance, in the context of the internship of Olivier Feudjio (Univ. de Paris), co-advised by Thibault Etienne and myself, we are currently applying our approach to new models and data sets, in order to study the role played by the localization of RNase E in mRNA stability.

One limitation of our approach, however, is the use of model simplification to remove the coupling of the ODEs introduced by the competition term. The price to pay is that the simplified model does no longer allow the analysis of the impact of extreme initial concentrations or Km values on bacterial mRNA decay. Although our study in [Etienne et al., 2020] suggests that this is not a key factor for the decay of most mRNAs, it would be interesting to study this specific point using the transcriptomic data. This requires additional developments of our approach that will be discussed in Chapter 8.

<footnote n="1">¹https://www.inrae.fr/actualites/comprendre-comment-bacteries-sadaptent-leur-environnement</footnote>
<figure><img src="image_39.png" /></figure>
CHAPTER 8

# Outlook

In this manuscript I have provided an overview of my past and current research activities at the interface of systems biology, bioinformatics, and microbiology. The material has been organized into five chapters around representative publications that I co-authored. As discussed in each chapter, the work has been carried out in the context of several research projects, in collaboration with a number of colleagues in France and abroad, and the co-supervision of PhD students and post-doctoral researchers. Additional research activities that have not been discussed in this manuscript relate to past and current collaborations with Andreas Kremling (TU München, Germany), Laurent Trilling (Univ. Grenoble Alpes), Tomas Gedeon (Montana State University, USA), Aline Métris and Jozsef Baranyi (formerly with Institute of Food Research, Norwich, UK), and Jean-Luc Gouzé (Inria Sophia Antipolis - Méditerranée). The evolution of the work described, from the qualitative to the quantitative modelling of biological systems, reflects the increased availability of quantitative data characterizing bacterial growth at multiple levels of cell organization. Most of my current and future research activities will be dedicated to the challenge of interpreting multiple heterogeneous data sets with mathematical models of bacterial growth and obtaining a more integrated view of cellular physiology. I describe these future directions of research below.

## 8.1 Genome-scale analysis of microbial physiology

### 8.1.1 Genome-scale analysis of cell metabolism

Foundations for this line of research have been laid by the on-going ANR project RIBECO (2018-2022), in which I am principal investigator for Inria Grenoble - Rhône-Alpes. In this project, we seek to learn more about the energetic burden imposed by the mRNA life cycle. Indeed, the continuing cycle of mRNA synthesis and degradation raises energetic constraints detrimental to cellular growth, in particular when the substrate consists of poor carbon sources. Our aim is to elucidate the connection between central carbon metabolism and RNA metabolism. While most of the work described in Chapter 6 relied on metabolomics data only, more data will be obtained in this project, such as transcriptomic and mRNA half-live data at the genome-wide level. This opens new opportunities to develop specific genome-scale models including both gene expression and metabolic pathways [O'Brien et al., 2013, Salvy and Hatzimanikatis, 2020]. Integration of these data sets is not an easy task as each comes with its own noise and bias, as seen in Chapter 4 for reporter gene and metabolomics data. I will therefore continue working on the problem of estimating biological quantities from primary data in collaboration with Eugenio Cinquemani and Hidde de Jong, as well as their implementation in user-friendly software.

Understanding the reprogramming of cell metabolism following a genetic or environmental perturbation is complicated. The main difficulty is to causally relate the predicted changes of fluxes to changes in gene expression and growth rate. A flux can vary directly with, for instance, the concentration of the enzyme catalysing the reaction and indirectly in response to the perturbation of another metabolic flux propagating through the network. I will therefore seek to develop approaches allowing one to identify at the genome-scale level the specific parts of metabolism directly affected by a perturbation. I have started investigating this question in collaboration with Marie-France Sagot and her team at Inria Lyon in the framework of the project MuSE, funded by the Complex
Systems Institute in Rhône-Alpes (IXXI), which I coordinated (2018-2020). We adapted a recently developed mixed-integer linear programming approach [Pusa et al., 2020], which we plan to generalize to the integration of multiple data sets (microarray and proteomics data, growth rates) and to apply to the problems studied in the context of RIBECO. Knowledge on RNA and energy metabolism obtained in this way should help us to design strategies that alleviate the energetic burden of the mRNA life cycle for biotechnological purposes. In particular, we will reengineer strains to help the efficient degradation of carbon sources derived from the pretreatment of agricultural and forestry residues.

The know-how developed for the analysis of metabolism and bacterial growth of *E. coli* can be applied to other organisms as well. I will work on two related projects in the coming months. For instance, Hidde de Jong and I have started a collaboration with Luiz de Carvalho from the Francis Crick Institute in London to analyse metabolic networks for the genus *Mycobacterium*, which we will reconstruct, and use these models to propose novel hypotheses on metabolic bottlenecks in the growth of mycobacteria. Another nascent project concerns the analysis of metabolic alterations during the development of the Parkinson disease, in collaboration with Florence Fauvelle from the Grenoble Institute of Neurosciences.

### 8.1.2 Genome-scale analysis of mRNA decay

The NLME framework generally yields good estimation results [Gonzalez et al., 2013] and it was critical in our integrative study of mRNA degradation. However, we had to overcome a number of issues relative to the noise of the data and identifiability problems, which may impede the extension of this approach to other high-throughput data sets taking into account additional biological processes, which is our objective in the RIBECO project. For instance in Equation 7.4, each reaction rate $v_i$ for a given mRNA *i* is a saturating function of the concentration and kinetic parameters of mRNA *i*, but also of the other mRNAs. This introduces a coupling between the ODEs resulting in correlations between model parameters, which could be intensified with the inclusion of other biological mechanisms. To circumvent the problem in our study, we simplified the model to decouple the equations. In collaboration with Aline Marguet, recently recruited in our group, and Eugenio Cinquemani, I will consider alternative NLME approaches by explicitly modelling the source of correlations introduced by the competition between mRNAs, based on an approach described in [Marguet et al., 2019]. Alternative approximations of the coupled models are also currently investigated in the context of the post-doctorate of Thibault Etienne.

Extending our approach to other data sets is not trivial, as each comes with its own noise. Another line of research will thus concern the development of appropriate preprocessing steps and error models in order to use heterogeneous data sets. Altogether this work should allow the exploitation of various time-series data sets of *E. coli* growth, metabolism, and gene expression obtained at TBI in the framework of the ANR RIBECO. A deeper understanding of the observed connections between mRNA degradation and the central carbon metabolism, and the regulatory mechanisms involved, would help us reach our final goal in the RIBECO project. That is, to propose changes in the life cycle of specific mRNAs and thus improve the degradation of vegetal biomass by microorganisms.

## 8.2 Resource allocation strategies in natural and engineered microorganisms

This direction of research is a development of our work in Chapter 5, in which we described a growth switch enabling the reallocation of nutrient resources from bacterial growth to the production of compounds. In order to understand the functioning of the growth switch on the molecular level, I am currently developing a mechanistic model of the gene expression machinery in *E. coli*, in collaboration with Hidde de Jong, Hans Geiselmann, and Jean-Luc Gouzé at Inria Sophia Antipolis - Méditerranée. In line with previous efforts in this direction [Dourado and Lercher, 2020, Weiße et al., 2015], the idea
is to provide a coarse-grained picture of the different macromolecular components, completed with the addition of RNA polymerase and its external control.

While many experimental studies monitored the macromolecular composition of *E. coli* cells growing at various steady-state growth rates, dynamical data for the gene expression machinery are scarce, if not absent. I will conduct experiments in the laboratory of Hans Geiselmann at LIPhy to provide us with original data on the dynamic adaptation of macromolecular cell composition to improve model estimation and prediction. A more detailed understanding of the growth switch will leverage the development of optimal strategies for producing metabolites or (heterologous) proteins, but on the fundamental level it may also provide novel insights into the dynamics of the adaptation of gene expression to environmental perturbations. Moreover, in the context of the RIBECO project, we expect the model to be useful for the design of strategies that tune the life cycle of specific mRNAs.

Another related line of research concerns the analysis of resource allocation strategies in bacteria. The foundations for this research axis have been laid by the ANR project Maximic (2017-2022), which involves modellers and experimentalists from our project-team, as well as other members of LIPhy, and specialists in control theory from Inria Sophia Antipolis - Méditerrannée (BIOCORE and McTAO project-teams), as well as Tomas Gedeon from Montana State University (USA). In this study, cells are viewed as self-replicators that try to grow optimally. The models in this case do not piece together all known biochemical reactions, but provide a coarse-grained picture of key cellular functions that captures the major fluxes of material and energy passing through the cell and fuelling growth. Such models may be instrumental for explaining a fundamental trade-off between rate and yield in the growth of microorganisms, that is, the fact that in microorganisms rapid growth generally comes at the cost of less efficient growth [Lipson, 2015].

## 8.3 From project-team IBIS to MICROCOSME

Inria project-teams have a maximal duration of twelve years. Our current project-team IBIS is hence terminating its life. With former permanent Ibis members (Hidde de Jong, Eugenio Cinquemani, Aline Marguet, Hans Geiselmann) and Muriel Cocaign-Bousquet, we are creating a new Inria project-team, MICROCOSME. I will take the lead of this new team, whose creation is currently under instruction. The on-going and future directions of research described are part of two of the four research axes of MICROCOSME. The other two axes are concerned with the analysis of the variability of bacterial growth, and heterogeneity within communities consisting of different microbial species and control of communities for biotechnological applications.

The start of a new scientific adventure...
<figure><img src="image_40.png" /></figure>
# Appendix
<figure><img src="image_41.png" /></figure>
# Curriculum Vitae

Delphine ROPERS

Born 29/11/1975 in Essey-lès-Nancy (54)

French nationality

Married, two children

## Current position

**Position:** Research scientist at Inria

**Professional address:** Centre de recherche Grenoble - Rhône-Alpes, Inovallée, 655 Avenue de l'Europe - CS 90051, 38334 Montbonnot Cedex, France

**Phone:** +33 4 76 61 53 72

**E-mail:** delphine.ropers@inria.fr

**Web:** https://team.inria.fr/ibis/delphine-ropers/

## Professional Experience

<table><tr><td>2017 -</td><td>Research scientist (Chargée de recherche classe normale) at Inria, Grenoble - Rhône-Alpes research centre (project-team Ibis)</td></tr><tr><td>2008 - 2017</td><td>Research scientist - first grade (Chargée de recherche 1<sup>re</sup> classe) at Inria, Grenoble - Rhône-Alpes research centre (project-team Ibis)</td></tr><tr><td>2006 - 2008</td><td>Research scientist - second grade (Chargée de recherche 2<sup>e</sup> classe) at Inria, Grenoble - Rhône-Alpes research centre (project-team Helix)</td></tr><tr><td>2003 - 2006</td><td>Post-doctoral researcher at Inria, Grenoble - Rhône-Alpes research centre (project-team Helix)</td></tr><tr><td>2002 - 2003</td><td>Doctoral researcher at CNRS Nancy</td></tr><tr><td>2001 - 2002</td><td>Temporary lecturer (ATER) at Université de Nancy I (96h)</td></tr><tr><td>1998 - 2001</td><td>Doctoral researcher at Université Nancy I</td></tr><tr><td>1998 - 2001</td><td>Instructor (monitorat) at Université Nancy I</td></tr></table>

## Academic Education

<table><tr><td>1998 - 2003</td><td>Doctoral thesis in Molecular and Cellular Biology at Laboratoire de Maturation des ARN et Enzymologie Moléculaire at CNRS/Université Nancy I. Title of the thesis: Experimental study of the role of SR proteins in the regulation of the HIV-1 virus RNA splicing, responsible of the human immunodeficiency, and mathematical modelling of these regulations.</td></tr><tr><td>1996 - 1998</td><td>Master of Science in Biochemistry at University of Nancy I. Title of the thesis: "Etude des éléments activant ou inhibant en cis l'épissage de l'ARN du virus VIH-1".</td></tr><tr><td>1993 - 1996</td><td>Bachelor of Science in Biochemistry at University of Nancy I.</td></tr></table>

## Scientific Production

### Articles in preparation

1. T.A. Etienne, L. Girbal, E. Cinquemani, M. Cocaign-Bousquet, **D. Ropers**. Integrative analysis of the physiological control of bacterial mRNA decay.
## Submitted articles

1. **D. Ropers**, Y. Couté, L. Faure, S. Ferré, D. Labourdette, A. Shabani, L. Trouilh, P. Vasseur, G. Corre, M. Ferro, M.-A. Teste, J. Geiselmann, H. de Jong. A multi-omics study of bacterial growth arrest in a synthetic biology application.

2. C. Roux, T.A. Etienne, E. Hajnsdorf, **D. Ropers**, A.J. Carpousis, M. Cocaïgn-Bousquet, L. Girbal. The essential role of mRNA degradation in understanding and engineering *E. coli* metabolism.

## Peer-reviewed articles

1. T.A. Etienne, M. Cocaïgn-Bousquet, **D. Ropers** (2020). Competitive effects in bacterial mRNA decay. *Journal of Theoretical Biology*, 504:110333.

2. M. Morin, B. Enjalbert, **D. Ropers**, L. Girbal, and M. Cocaïgn-Bousquet (2020). Genome-wide stabilization of mRNA during a 'feast-to-famine' growth transition in *Escherichia coli*. *mSphere*, 5:e00276-20.

3. S. Pinhal, **D. Ropers**, J. Geiselmann, H. de Jong. Acetate metabolism and the inhibition of bacterial growth by acetate (2019). *Journal of Bacteriology*, 201:e00147-19.

4. I. Belgacem, S. Casagranda, E. Grac, **D. Ropers**, J.-L. Gouzé (2018). Reduction and stability analysis of a transcription-translation model of RNA polymerase. *Bulletin of Mathematical Biology*, 80(2), 294-318.

5. A. Kremling, J. Geiselmann, **D. Ropers**, H. de Jong (2018). An ensemble of mathematical models showing diauxic growth behaviour. *BMC Systems Biology*, 12:82.

6. S. Casagranda, S. Touzeau, **D. Ropers**, J.-L. Gouzé (2018). Principal process analysis of biological models. *BMC Systems Biology*, 12(1):68.

7. H. de Jong, S. Casagranda, N. Giordano, E. Cinquemani, **D. Ropers**, J. Geiselmann, J.-L. Gouzé (2017). Mathematical modeling of microbes: Metabolism, gene expression, and growth. *Journal of the Royal Society Interface*, 14:20170502.

8. M. Morin, **D. Ropers**, E. Cinquemani, J.C. Portais, B. Enjalbert, M. Cocaïgn-Bousquet (2017). The Csr system regulates *Escherichia coli* fitness by controlling glycogen accumulation and energy levels. *mBio*, 8(5): e01628-17.

9. H. de Jong, **D. Ropers**, J. Geiselmann (2017). Resource reallocation in bacteria by reengineering the gene expression machinery. *Trends in Microbiology*, 25(6):480-493.

10. E. Cinquemani, V. Laroute, M. Cocaïgn-Bousquet, H. de Jong, D. Ropers (2017). Estimation of time-varying growth, uptake and excretion rates from dynamic metabolomics data. *Bioinformatics*, 33(14):i301-i310.

11. A. Métris, S.M. George, **D. Ropers** (2017). Piecewise linear approximations to model the dynamics of adaptation to osmotic stress by food-borne pathogens. *International Journal of Food Microbiology*, 240:63-74.

12. **D. Ropers** and A. Métris. Osmotic stress response to NaCl in *Escherichia coli*: qualitative modeling and simulation data (2016), *Data in Brief*, 9: 606-612.

13. M. Morin, **D. Ropers**, F. Létisse, S. Laguerre, J.C. Portais, M. Cocaïgn-Bousquet, B. Enjalbert (2016). The post-transcriptional regulatory system CSR controls the balance of metabolic pools in upper glycolysis of *Escherichia coli*. *Molecular Microbiology*, 100(4):686-700.

14. J. Izard, CDC Gomez Balderas, **D. Ropers**, S. Lacour, X. Song, Y. Yang, AB Lindner, J. Geiselmann, H. de Jong (2015). A synthetic growth switch based on controlled expression of RNA polymerase *Molecular Systems Biology*, 11(840).
15. A. Kremling, J. Geiselmann, **D. Ropers**, H. de Jong, Understanding carbon catabolite repression in *Escherichia coli* using quantitative models (2015). *Trends in Microbiology*, 23(2):99-109.

16. V. Zulkower, M. Page, **D. Ropers**, J. Geiselmann, H. de Jong (2015), Robust reconstruction of gene expression profiles from reporter gene data using linear inversion. *Bioinformatics*, 15:31(12):i71-9

17. M. Trauchessec, M. Jaquinod, A. Bonvalot, V. Brun, C. Bruley, **D. Ropers**, H. de Jong, J. Garin, G. Bestel-Corre, M. Ferro (2014), Mass spectrometry-based workflow for accurate quantification of *E. coli* enzymes : how proteomics can play a key role in metabolic engineering, *Molecular and Cellular Proteomics*. 13(4):954-968.

18. G. Baptist, C. Pinel, C. Ranquet, J. Izard, **D. Ropers**, H. de Jong, J. Geiselmann (2013), A genome-wide screen for identifying all regulators of a target gene, *Nucleic Acids Research*, 41(17):e164.

19. S. Berthoumieux, H. de Jong, G. Baptist, C. Pinel, C. Ranquet, **D. Ropers**, J. Geiselmann (2013), Shared control of gene expression in bacteria by transcription factors and global physiology of the cell, *Molecular Systems Biology*, 9:634. Editors' choice in Science

20. V. Baldazzi, **D. Ropers**, J. Geiselmann, D. Kahn, H. de Jong (2012), Importance of metabolic coupling for the dynamics of gene expression following a diauxic shift in *Escherichia coli*, *Journal of Theoretical Biology*, 295: 100-115.

21. G. Batt, B. Besson, P.-E. Ciron, H. de Jong, E. Dumas, J. Geiselmann, R. Monte, P.T. Monteiro, M. Page, F. Rechenmann, **D. Ropers** (2012), Genetic Network Analyzer: A tool for the qualitative modeling and simulation of bacterial regulatory networks, J. van Helden, A. Toussaint, D. Thieffry (eds), Bacterial Molecular Networks: Methods and Protocols, *Methods in Molecular Biology*, Humana Press, Springer, New York, 439-462.

22. P.T. Monteiro, P.J. Dias, **D. Ropers**, A.L. Oliveira, I. Sá-Correia, M.C. Teixeira and A.T. Freitas (2011), Qualitative modeling and formal verification of the FLR1 gene mancozeb response in *Saccharomyces cerevisiae*, *IET Systems Biology*, 5(5): 308-316.

23. **D. Ropers**, V. Baldazzi, H. de Jong (2011), Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in *Escherichia coli*, *IEEE/ACM Transactions on Computational Biology and Bioinformatics*, 8(1):166-181.

24. V. Baldazzi, **D. Ropers**, Y. Markowicz, D. Kahn, J. Geiselmann, H. de Jong (2010), The carbon assimilation network in *Escherichia coli* is densely connected and largely sign-determined by directions of metabolic fluxes, *PLoS Computational Biology*, 6(6):e1000812.

25. H. de Jong, C. Ranquet, **D. Ropers**, C. Pinel, J. Geiselmann (2010), Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria, *BMC Systems Biology*, 4:55.

26. F. Boyer, B. Besson, G. Baptist, J. Izard, C. Pinel, **D. Ropers**, J. Geiselmann, H. de Jong (2010), A MATLAB program for the analysis of fluorescence and luminescence reporter gene data, *Bioinformatics*, 26(9):1262-1263.

27. F. Corblin, S. Tripodi, E. Fanchon, **D. Ropers**, L. Trilling (2009), A declarative constraint-based method for analyzing discrete genetic regulatory networks, *BioSystems*, 98(2):91-104.

28. J.-M. Saliou, C.F. Bourgeois, L. Ayadi-Ben Mena, **D. Ropers**, S. Jacquenet, V. Marchand, J. Stévenin and C. Branlant (2009), Role of RNA structure and protein factors in the control of HIV-1 splicing, *Frontiers in Biosciences*, 14:2714-2729.