Resumé en français

Cette thèse de doctorat porte sur l'analyse du métabolisme des micro-organismes. Le métabolisme est l'ensemble des réactions d'un organisme. Il peut être modélisé comme un réseau métabolique qui contient les métabolites présents dans un organisme et les réactions correspondantes qui les transforment. Les réseaux métaboliques peuvent être, par exemple, représentés sous forme de hypergraphes ou de matrices stoechiométriques. La modélisation du métabolisme peut être utilisée pour obtenir des informations sur l'activité métabolique d'un micro-organisme et pour prédire certains comportements.

La disponibilité croissante des données métabolomiques et leur analyse améliorent la compréhension des mécanismes cellulaires. L'intégration des données métabolomiques dans les réseaux métaboliques permet de comprendre comment les systèmes biologiques réagissent à différentes perturbations. Cela peut être, par exemple, l'adaptation du micro-organisme à un changement de milieu, au stress ou l'impact des perturbations génétiques sur l'activité métabolique du micro-organisme. La modification de l'activité métabolique d'un micro-organisme est appelée changement métabolique (en anglais "metabolic shift"). Comme il est important de comprendre les changements métaboliques dans différentes conditions pour faire progresser la recherche dans différents domaines comme la bioingénierie ou la santé humaine, il est nécessaire de disposer de méthodes de calcul qui facilitent la compréhension des données métabolomiques disponibles.

Un autre sujet commun dans l'analyse des réseaux métaboliques est le calcul des stratégies d'intervention optimales. Les micro-organismes sont déjà utilisés à l'échelle industrielle pour produire des substances chimiques cibles importantes. Il est essentiel de comprendre comment les différentes parties du métabolisme interagissent et sont liées entre elles pour trouver des moyens d'améliorer les taux de production. Le métabolite souhaité peut souvent n'être qu'un sous-produit du métabolisme habituel lorsque le micro-organisme est en croissance et son rendement peut être très faible. Il est donc nécessaire de modifier le métabolisme du micro-organisme afin d'augmenter le rendement du composé cible souhaité et établir ainsi une production plus efficace. L'un des domaines de la biologie synthétique est le génie génétique des micro-organismes pour améliorer la production de substances chimiques cibles importantes par un micro-organisme. L'idée est d'éliminer des réactions spécifiques en manipulant les gènes métaboliques qui codent pour les enzymes catalytiques. En fonction du rôle que les réactions d'élimination jouent dans le métabolisme, l'organisme peut utiliser différentes voies pour compenser son manque. Les réactions éliminées doivent être choisies de manière que le métabolisme modifié entraîne une augmentation de la production de la substance chimique cible.

En raison de la structure complexe du métabolisme, il est nécessaire de modéliser les perturbations in silico afin de proposer des stratégies d'intervention qui pourraient conduire au meilleur résultat in vivo. En prédisant et en analysant le comportement modifié et en identifiant les meilleures perturbations à l'aide d'approches computationnelles, la conception d'expériences pratiques peut être guidée.

De plus en plus d'approches sont mises au point pour calculer des knock-outs optimaux pour des scénarios différents. Certaines ne sont applicables qu'à des modèles de réseaux plus petits qui ne représentent qu'une version condensée ou que certaines parties du métabolisme. L'application aux réseaux métaboliques à l'échelle du génome qui modélisent le métabolisme complet d'un microorganisme de manière aussi détaillée que possible reste particulièrement difficile. En raison de leur taille et de leur complexité, l'effort de calcul des méthodes appliquées peut augmenter considérablement.

Il y a cependant un aspect qui n'est souvent pas pris en compte lors du calcul des knock-outs optimaux qui devraient augmenter la production cible. Dans certains cas, le produit chimique cible est un sous-produit qui est en fait toxique pour le micro-organisme. Par conséquent, l'accumulation du métabolite cible peut inhiber la croissance du micro-organisme. Elle peut également diminuer le taux de production, ce qui rendra l'application à l'échelle industrielle moins efficace.

Les knock-outs restreignent généralement le métabolisme d'un micro-organisme parce qu'ils enlèvent certaines fonctionnalités. Par conséquent, dans le cadre mentionné ci-dessus, il est crucial de s'assurer que les knock-outs effectués n'inhibent pas des processus importants que le micro-organisme peut utiliser pour établir une tolérance contre la cible toxique. Cela ne contribuera pas spécifiquement à augmenter la résistance contre le produit chimique cible, mais cela permet de s'assurer que la tolérance naturelle est préservée.

J'ai déjà travaillé sur les réseaux métaboliques pendant mon mastère, ce qui a également été l'une des raisons pour lesquelles je me suis intéressée à ce doctorat. J'ai travaillé sur deux approches pour les réseaux métaboliques qui utilisent l'optimisation multi-objectifs. Une approche était applicable aux communautés microbiennes, c'est-à-dire pour un scénario où différents micro-organismes interagissent les uns avec les autres. La seconde utilisait l'optimisation multi-objectifs pour identifier les compromis possibles entre la production de biomasse et la production cible dans les micro-organismes. Des knock-outs sont énumérées consécutivement. Cette deuxième approche a donné une première idée de base pour l'approche développée qui sera la partie principale de cette thèse. Elle sera présentée dans le chapitre 2.

Au début de cette thèse, dans le chapitre 1, je ferai un bref résumé sur certains sujets qui sont importants pour la compréhension des deux approches qui seront présentées dans les chapitres 2 et 3. Je ferai d'abord une brève introduction à la programmation linéaire, en particulier à la programmation linéaire multi-objectifs qui sera utilisée dans l'approche décrite au chapitre 2. Ensuite, je donnerai une introduction aux réseaux métaboliques, à leurs représentations et aux méthodes associées. Dans les deux dernières sections, je résumerai les approches communes pour prévoir les stratégies d'intervention et analyser les changements métaboliques.

Le chapitre 2 représente le principal travail que j'ai effectué pendant mon doctorat. Cette partie se concentre sur une approche d'exploitation du métabolisme des micro-organismes. La motivation principale était de développer une stratégie pour prédire les knock-outs optimaux qui augmentent la production d'une substance chimique cible dans un micro-organisme dans le cas où le métabolite cible produit est toxique pour le micro-organisme utilisé. En outre, je souhaitais une approche qui soit également applicable aux réseaux métaboliques à l'échelle du génome.

L'approche développée se compose de deux parties différentes. Dans la première partie, un problème d'optimisation multi-objectifs est formulé qui calcule les compromis entre la production de biomasse, la production cible et un score qui mesure la résistance possible à la toxicité contre le métabolite cible toxique. Dans la deuxième partie, on énumère des knock-outs qui devraient imposer des distributions de flux spécifiques. Dans un premier temps, un MILP a été proposé pour énumérer les différents knock-outs. En raison de ses limites et de ses défauts, il était nécessaire de développer une deuxième idée qui est basée sur l'identification et l'isolation de sous-réseaux.

La méthode proposée est appliquée à l'étude de cas de la production d'éthanol dans la levure. L'éthanol est déjà produit par la levure dans l'industrie et il y a un intérêt croissant pour l'éthanol en raison de son utilisation comme biocarburant. Cependant, un facteur limitant important pour la production d'éthanol dans la levure est en effet la toxicité de l'éthanol pour la levure. Lorsque l'éthanol s'accumule, il inhibe la croissance mais entraîne également un déclin de la production d'éthanol.

Il pourrait donc être important de tenir compte de cet aspect lorsque l'on essaie d'améliorer la production d'éthanol en introduisant des knock-outs qui entraîneront un changement dans l'activité métabolique de la levure. En appliquant notre approche à l'étude de cas de la production d'éthanol dans la levure, nous avons pu calculer des ensembles de désactivation avec moins de 20 réactions. Il reste à déterminer quelle est la plus petite taille possible. En outre, nous avons besoin d'une évaluation biologique plus poussée pour les ensembles de désactivation proposés afin d'établir la viabilité des réactions suggérées dans la pratique. Les avantages de notre approche sont qu'elle est applicable sur des réseaux métaboliques à l'échelle du génome comme nous l'avons montré en utilisant le modèle de la levure 5.01. En outre, le cadre est flexible et il devrait être possible de modifier les objectifs pour différents scénarios. Nous sommes donc également intéressés par l'application de notre approche à d'autres exemples afin de confirmer son adaptabilité. Nous prévoyons de peaufiner l'approche et de la rendre disponible comme outil sur le Gitlab de l'équipe. Après avoir obtenu un avis biologique supplémentaire des collaborateurs sur nos résultats, nous voulons soumettre ce travail sous forme d'article avant la fin de l'année. Dans le chapitre 3, je présenterai une nouvelle méthode de calcul qui s'est concentrée sur étudier l'activité métabolique des micro-organismes. Cette approche est appelée Totoro, ce qui est l'abréviation de "Transient respOnse to meTabOlic pertuRbation inferred at the whole netwOrk level". Cette approche est le résultat d'une collaboration avec Ricardo Andrade et Mariana Ferrarini. Un ancien membre du groupe, Alice Julien-Laferrière, a déjà travaillé sur cette approche lors de son doctorat. Nous avons soumis un papier sur ce travail à Bioinformatique en Septembre 2020.

Totoro intègre les concentrations internes de métabolites qui ont été mesurées avant et après une perturbation dans des reconstructions métaboliques à l'échelle du génome. Il prédit les réactions qui étaient actives pendant l'état transitoire qui a suivi la perturbation. Il s'agit d'une approche basée sur les contraintes qui prend en compte la stoechiométrie du réseau. Elle minimise le changement des concentrations pour les métabolites non mesurés et également le nombre de réactions actives pendant l'état transitoire pour tenir compte d'une hypothèse parcimonieuse. Totoro est un outil disponible librement. Il a été implémenté en C++ et peut être consulté à l'adresse https://gitlab.inria. fr/erable/totoro. Nous avons appliqué notre méthode à trois expériences d'impulsions dans Escherichia coli pour montrer qu'elle peut récupérer des voies actives connectées et prédire des directions distinctes pour des réactions réversibles qui sont conformes aux observations biologiques connues. Nous avons utilisé un modèle de base et un modèle à l'échelle du génome de E. coli pour montrer que notre approche est également applicable à des modèles de réseaux d'une taille plus grande. Un autre projet sur lequel j'ai commencé à travailler pendant mon doctorat est une collaboration avec Marianne Borderes. Nous développons une approche visant à combiner les résultats de plusieurs méthodes de regroupement appliquées aux données de la métagénomique. Ce projet ne sera pas présenté dans cette thèse car il est toujours en cours et, de plus, il n'est pas directement lié à l'analyse des réseaux métaboliques. This PhD thesis is about the analysis of the metabolism of microorganisms. The metabolism is the set of reactions of an organism. It can be modeled as a metabolic network which contains the metabolites that are present in an organism and the corresponding reactions that transform them. Metabolic networks can be, for instance, represented as a hypergraph or a stoichiometric matrix.

Modeling the metabolism can be used to gain insights on the metabolic activity of a microorganism and to predict certain behaviors,

The increasing availability of metabolomic data and their analysis are improving the understanding of cellular mechanisms. Integrating metabolomic data into metabolic networks makes it possible to understand how biological systems respond to different perturbations. This can be, for example, the adaption of the microorganism to a change in the medium or to stress, or else the impact of genetic perturbations on its metabolic activity. The change in the metabolic activity of a microorganism is called metabolic shift. Since understanding metabolic changes under different conditions is important for advancing research in different fields like bioengineering or human health, there is a need for computational methods that facilitate the comprehension of available metabolomic data.

Another common topic in the analysis of metabolic networks is the computation of optimal intervention strategies. Microorganisms are already used on an industrial scale to produce important target chemicals. Understanding how different parts of the metabolism interact and are linked together is crucial to find ways that will improve production rates. The desired compound might often just be a by-product of the usual metabolism when the microorganism is growing and its yield can be very low. Thus, it is necessary to modify the metabolism of the microorganism to increase the yield of the desired target compound and to establish a more efficient production. One of the areas of synthetic biology is the genetic engineering of microorganisms to improve the production of important target chemicals by a microorganism. The idea is to knockout specific reactions by manipulating the metabolic genes that code for the catalyzing enzymes. Depending on the role the knocked out reactions play in the metabolism, the organism might use different pathways to compensate for its lack. The reactions that are knocked out must be chosen in a way that the altered metabolism leads to an increase in the production of the target chemical.

Due to the complex structure of metabolisms, there is a need to model perturbations in silico in order to propose intervention strategies that might lead to the best outcome in vivo. By predicting and analyzing the altered behavior and identifying optimal knockouts with computational approaches, the design of practical experiments can be supported and guided.

More and more approaches are developed that aim to compute optimal knockouts for different scenarios. Some are only applicable to smaller network models which are only representing a condensed version or certain parts of the metabolism. Especially the application to genome-scale metabolic networks which model the whole metabolism of a microorganism as detailed as possible remains challenging. Due to their size and complexity, the computational effort of the applied methods can increase significantly.

When computing optimal knockouts that should increase the target production, there is however one aspect that is often not taken into account. In some cases, the desired target chemical is a by-product that is actually toxic for the microorganism and, therefore, its accumulation can severely inhibit the growth of the microorganism. It can also decrease the production rate which will make the application on a industrial scale less efficient.

Knockouts usually restrict the metabolism of a microorganism because they are taking away certain functionalities. Hence, in the above-mentioned setting, it is crucial to ensure that the inserted knockouts do not inhibit important processes that a microorganism can use to establish a tolerance against the toxic target. This will not specifically help to increase the resistance against the target chemical but it ensures that the natural tolerance is preserved.

I already worked on metabolic networks during my master's degree which was also one of the reasons why I got interested in this PhD. I worked on two approaches for metabolic networks that use multiobjective optimization. One approach was applicable to microbial communities, i.e., for a scenario where different microorganisms interact with each other. The second one was using multi-objective optimization to identify possible tradeoffs between biomass production and target production in microorganisms. Knockouts for each tradeoff were enumerated afterwards. This second approach gave a first, basic idea for the developed approach that will be the main part of this thesis. It will be presented in Chapter 2.

In the beginning of this thesis, in Chapter 1, I will give a brief summary on certain topics that are important for the understanding of the two approaches that will be presented in the Chapters 2 and 3. I will first do a short introduction to linear programming, especially multi-objective linear programming which will be used in the approach described in Chapter 2. Afterwards, I will give an introduction to metabolic networks, their representations and associated methods. In the last two sections, I will summarize common approaches to predict intervention strategies and to analyze metabolic shifts.

Chapter 2 represents the main work that I did during my PhD. This part focuses on an approach to exploit the metabolism of microorganisms. The main motivation was to developed a strategy to predict optimal knockouts that increase the production of a target chemical in a microorganism in the scenario where the produced target metabolite is toxic for the utilized microorganism. Furthermore, I aimed for an approach that is also applicable to genome-scale metabolic networks.

The developed approach consists of two different parts. In the first part, a multi-objective optimization problem is formulated that computes tradeoffs between biomass production, target production and a score that measures the possible toxicity resistance against the toxic target metabolite. In the second part, knockouts are enumerated that should enforce specific flux distributions. As a first idea, a mixed-integer linear program was proposed to enumerate different knockouts. Due to its limitations and flaws, there was a need to develop a second idea which is based on identifying and cutting off subnetworks.

The proposed method is applied to the case-study of ethanol production in yeast. Ethanol is already produced by yeast in industry and there is a growing interest in ethanol due to its use as bio-fuel. However, an important limiting factor for the production of ethanol in yeast is indeed the toxicity of ethanol for yeast. When ethanol accumulates, it inhibits growth but leads also to a decline of the ethanol production. It might thus be important to consider this aspect when trying to improve the ethanol output by introducing knockouts that will lead to a change in the metabolic activity of yeast.

Based on this case-study, I will describe insights that we could gain from both parts of the approach, especially on the identified subnetworks. We plan on submitting a paper on this subject this year after obtaining some more biological evaluation of the obtained results. Furthermore, the developed approach will be made available as tool on the Gitlab of our group.

In Chapter 3, I will present a novel computational method that focused on investigating the metabolic activity of microorganisms. The approach is called Totoro which is short for "Transient respOnse to meTabOlic pertuRbation inferred at the whole netwOrk level". This approach is the result of a collaboration with Ricardo Andrade and Mariana Ferrarini. A former member of the group, Alice Julien-Laferrière, already worked on this approach during her PhD. We submitted a paper on this work to Bioinformatics in September 2020.

Totoro integrates internal metabolite concentrations that were measured before and after a perturbation into genome-scale metabolic reconstructions. It predicts reactions that were active during the transient state that occurred after the perturbation. It is a constraint-based approach that takes the stoichiometry of the network into account. It minimizes the change in concentrations for unmeasured metabolites and also the number of active reactions during the transient state to account for a parsimonious assumption. Totoro is a freely available tool implemented in C++ and it can be accessed at https://gitlab.inria.fr/erable/totoro. We applied our method to three published pulse experiments in Escherichia coli to show that it can retrieve connected active pathways and predict distinct directions for reversible reactions that are in accordance with known biological observations. We used a core model and a genome-scale model of E. coli to further demonstrate that our approach is also applicable to larger network models.

Another project that I started working on during my PhD is a collaboration with Marianne Borderes.

We are developing an approach to combine the results of multiple clustering methods applied to metagenomics data. This project will not be presented in this thesis because it is still ongoing and additionally, it is not directly connected to the analysis of metabolic networks.

1 Theoretical background In this chapter, I will first present a short introduction to linear programming, especially multiobjective linear programming which will be used in the approach described in Chapter 2. Afterwards, I will give an introduction to metabolic networks, their representations and associated methods. In the last two sections, I will summarize common approaches to predict intervention strategies and to analyze metabolic shifts.

Linear programming

As described in [START_REF] Dantzig | Linear programming 1: introduction[END_REF], linear programming maximizes or minimizes a linear objective function for the variables that are subject to linear constraints. Introductions to linear programming and its applications are, for example, given in [START_REF] Bertsimas | Introduction to linear optimization[END_REF][START_REF] Dantzig | Linear programming and extensions[END_REF][START_REF] Schrijver | Theory of linear and integer programming[END_REF][START_REF] Gass | Linear programming: methods and applications[END_REF]. Linear programming is part of the larger field of mathematical programming that includes amongst others integer programming and nonlinear programming.

Problem formulation

Mathematically, a linear program finds the values for the problem variables x that maximize an objective function z (Equation 1. 1). The problem variables are subject to m linear constraints (Equation 1.2).

max x z = f (x) = c T x (1.1) s.t. Ax ≤ b (1.2) A ∈ R m,n , x ∈ R n , b ∈ R m , c ∈ R n (1.3)
A vector x ∈ R n is a feasible solution if it does not violate any of the given constraints. The feasible region is the set of all feasible solutions. Geometrically, the feasible region is described by a polyhedron. A solution x * is optimal if it it a feasible solution and f

(x * ) = max{f (x)|Ax ≤ b, x ∈ R n }.
Linear programs can be infeasible which means that no solution exists that satisfies all constraints. Furthermore, it can have several optimal solutions which means that there are multiple feasible solutions that lead to the same optimal value for the objective function. The linear program is unbounded if the value for the objective function can be increased (or decreased) without any limit. If all variables are integer, the resulting problem is called integer linear program (ILP). If some are integer and some are continuous, it is a mixed integer linear program (MILP). Solvers for linear programs are, for example, IBM ILOG CPLEX (IBM, 2019), Gurobi Optimization (Gurobi Optimization, 2020) or SCIP [START_REF] Achterberg | Constraint integer programming: A new approach to integrate cp and mip[END_REF][START_REF] Achterberg | Scip: solving constraint integer programs[END_REF].

Multi-objective linear programming

In contrast to a linear program which optimizes one objective function, a multi-objective linear program minimizes or maximizes k objective functions at the same time (Equation 1.4). The different objective functions can be contradictory. Multi-objective programs are, for instance, described in [START_REF] Ehrgott | Multicriteria Optimization[END_REF] and [START_REF] Antunes | Multiobjective linear and integer programming[END_REF].

max x z 1 = f 1 (x) . . . (1.4) max x z k = f k (x) s.t. Ax ≤ b (1.5) x ∈ R n (1.6)
A multi-objective linear program has, like a single-objective linear program, linear constraints that define the set of feasible solutions X in the decision space (Equation 1. [START_REF]coli core model -Results for glucose pulse (λ = 0.5, = 5)[END_REF]. In contrast to a linear program, here, the objective space is not one but k-dimensional where k matches with the number of objective functions. If k is one, the optimization problem is a single-objective linear program. The feasible set Z in the objective space is defined by all the points which are part of the feasible set in the decision space (Equation 1.8).

X = {x ∈ R n |Ax ≤ b, x j ∈ Z, j ∈ I} (1.7) Z = {z ∈ R k |z i = f i (x), x ∈ X, i = 1, ..., k} (1.8)
The concept of optimality in multi-objective optimization problems is different from single-objective linear programs because it is possible that the objective functions conflict with each other. This means that it is not necessarily possible that all objective functions can reach their optimal value at the same time. Solutions that are optimal in the sense of multi-objective programs are called efficient. Hence, a new relation is introduced to enable the comparison of two solutions z 1 , z 2 ∈ R k of a multi-objective linear program (Equation 1.9). This relation can be used to decide which solution is better or which solution dominates the other. A solution z 1 dominates another solution

z 2 if z 1 z 2 . z 1 z 2 ⇐⇒ z 1 i ≥ z 2 i , i = 1, ..., k and z 1 = z 2 (1.9)
A solution x * is called efficient if the corresponding point in the objective space is not dominated by any other point. The set X E includes all efficient solutions (Equation 1.10). The corresponding points in the objective space are nondominated points because there exit no other points that dominate them. The set Z N contains all nondominated points (Equation 1. [START_REF]coli core model -First ten solutions for glucose pulse (λ = 0.1[END_REF].

X E = {x ∈ X| x ∈ X : f (x) f (x)} (1.10) Z N = {z ∈ Z|z = f (x), x ∈ X E } (1.11)
In other words, a solution is treated as an efficient solution if it is not possible to increase the value of one objective function further without decreasing the values of one or more of the other objective functions. The difference between nondominated and dominated points is further illustrated in Figure 1.1a. The set of all nondominated points is called Pareto front which is why efficient or nondominated solutions are also referred to as Pareto optimal or Pareto efficient solutions. Furthermore, it is possible to distinguish between supported and unsupported nondominated solutions. Unsupported nondominated solutions are dominated by a (infeasible) convex combination of other nondominated points (Figure 1.1b).

Concepts to solve multi-objective linear programs

In [START_REF] Antunes | Multiobjective linear and integer programming[END_REF] and [START_REF] Ehrgott | Multicriteria Optimization[END_REF], some basic techniques to solve multi-objective linear programs are given. They include the weighted-sum method and the -constraint method which both rely on transforming the multi-objective linear program into several single-objective linear programs that have to be solved separately.

In the case of the weighted-sum method, the corresponding single-objective linear programs maximize the weighted sum for all objective functions of the multi-objective linear program (Equation 1.12. By changing the weights λ in the sum and resolving the changed single-objective linear program, different efficient solutions for the multi-objective linear program can be computed. Usually, the weights are normalized (Equation 1.13).

max x k i=1 λ i f i (x) (1.12) k i=1 λ i = 1 (1.13)
The weighted-sum method has the advantage that no additional constraints have to be added when transforming the problem. The computation complexity is therefore not changed and the transform problem requires the same computational effort as the single-objective version of the problem. On the other hand, it is only possible to compute supported nondominated points. Unsupported nondom-

z 1 z 2 a b c d e f 2 2 2 1 Figure 1.2: Basic concept of the -constraint method for k = 2.
The idea for this figure is taken from [START_REF] Antunes | Multiobjective linear and integer programming[END_REF]. In this example, z 1 is kept as objective function in the transformed problem. The remaining objective function z 2 is turned into the constraint z 2 ≥ 2 . In a first step, 2 can be set to zero, that is only z 1 is considered and maximized which means that the nondominated point f can be found. Afterwards, the value of 2 needs to be increased to compute other efficient solutions.

inated points cannot be found which means that in nonconvex problems, not all efficient solutions can be found by this approach. Similar to the weighted-sum method, the -constraint method also transforms the problem into multiple single-objective linear programs. In the transformed problem, only one of the objective functions is kept (Equation 1.14). The other objective functions are transformed into constraints (Equation 1.15).

max x f j (x) (1.14) s.t. f i (x) ≥ i , i = 1, ..., k, i = j (1.15)
Changing the bounds for i allows to compute different efficient solutions for the multi-objective linear program. An illustration of this idea can be found in Figure 1.2 for the case of two objective functions. In contrast to the weighted-sum method, the -constraint method can find all efficient solutions but as a result of the additional constraints, the transformed problem can be more difficult to solve. Another disadvantage that both problems have in common is that it is not clear when or if all nondominated points have been found.

PolySCIP [START_REF] Borndörfer | Polyscip[END_REF] is, for example, a solver for multi-objective linear programs.

Metabolic networks

The metabolism is the set of chemical reactions that are taking place in an organism. The metabolism can be modeled through a metabolic network. An introduction to metabolic networks, their definitions and also associated methods are given in [START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF][START_REF] Klamt | Stoichiometric and constraint-based analysis of biochemical reaction networks[END_REF][START_REF] Cottret | Graph methods for the investigation of metabolic networks in parasitology[END_REF]. The following definitions are taken from [START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF]. Metabolic networks consist of chemical compounds, biochemical reactions, enzymes and also genes.

The chemical compounds are also called metabolites. These are small molecules inside an organism that can be imported and exported but also synthesized and degraded. A biochemical reaction pro-duces a set of metabolites, the products, from another set of metabolites, the substrates. Reactions can import metabolites from an external source and they convert them, for example, into other building blocks needed by the organism to survive and to grow. Metabolites that are not needed by the organism and considered as waste can be excreted. Theoretically, reactions can take place in both directions which means that the set of products and substrates are interchangeable. These reactions are therefore reversible. However, depending on the physiological conditions, there are reactions that only take place in one specific direction. These reactions are called irreversible.

Reactions are catalyzed by enzymes. There are some reactions that happen spontaneously and that can be accelerated by enzymes but most reactions need to be catalyzed by enzymes or they cannot take place. Enzymes are proteins or protein complexes that are coded by one or by multiple genes.

Understanding the relationships between genes, enzymes and reactions in detail is difficult because a single enzyme can catalyze different reactions and a single reaction can be catalyzed by different enzymes. Moreover, sometimes co-factors are necessary to enable the catalysis of a reaction by enzymes. Co-factors are small molecules that can bind enzymes thereby changing the activity of the enzyme.

A metabolic network connects metabolites with reactions that transform them. Source reactions can uptake necessary metabolites from external sources and sink reactions excrete waste products or excess metabolites. Certain sets of reactions as a whole are referred to as metabolic pathway if they are part of the global synthesis or degradation of a specific metabolite with intermediate steps. For example, the glycolysis converts glucose into pyruvate and can be seen as one metabolic pathway.

Metabolic network reconstruction and databases

Metabolic networks are reconstructed based on the available knowledge on relations between genes, enzymes and reactions [START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF]. Usually, the reconstruction depends on comparative genomics but also on the use of metabolomic data which quantifies the metabolites that are present in the given organism. It is possible to infer from comparative genomics certain reactions that are present in an organism. In a first step, the functional annotation for genes is used to determine which enzymes are existent. Afterwards, the functional annotation is linked to the reaction by identifying which reactions are catalyzed by the present enzymes. Metabolic networks can be reconstructed automatically but the results can be erroneous. Due to their low quality, they should be used with care. To obtain accurate reconstructions, manual corrections based on the literature are necessary.

Metabolic pathways are available in databases like KEGG [START_REF] Kanehisa | KEGG: kyoto encyclopedia of genes and genomes[END_REF] and BioCyc [START_REF] Karp | The biocyc collection of microbial genomes and metabolic pathways[END_REF][START_REF] Caspi | The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases[END_REF]. Some databases are specialized on a specific organism, e.g.

EcoCyc [START_REF] Karp | The ecocyc database[END_REF] for Escherichia coli or HumanCyc [START_REF] Romero | Computational prediction of human metabolic pathways from the complete human genome[END_REF] for the human metabolism. Depending on the organism of interest, it can be difficult to find accurate metabolic networks, especially genome-scale metabolic networks which represent the complete metabolism of an organism. Problematically, the catalogued metabolic pathways can also differ between databases [START_REF] Altman | A systematic comparison of the MetaCyc and KEGG pathway databases[END_REF].

Graph representation

Metabolic networks can be represented in different ways. One of the most common representations is the modeling of a metabolic network as a graph [START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF]. In general, a graph G is formally defined as a pair (V, E). V corresponds to the set of vertices and E to the set of edges, where an edge is a subset of 2 vertices in V . The choice of which graph model should be used depends on the application. The compound graph models the metabolic network in a way where the compounds are represented by the vertices of the graph and the edges correspond to the reactions. If an edge connects two compounds, it means that there exists a reaction that transforms one compound into the other one. Therefore, one compound is a substrate and the other one is a product (see Figure 1.3a).

The problem is that there are many reactions that have multiple substrates and multiple products.

In order to represent these connections, hypergraphs can be used to describe the metabolic network structure [START_REF] Yeung | Estimation of the number of extreme pathways for metabolic networks[END_REF][START_REF] Pearcy | Hypergraph models of metabolism[END_REF][START_REF] Klamt | Hypergraphs and cellular networks[END_REF].

Hypergraphs are described in [START_REF] Berge | Graphs and Hypergraphs[END_REF]. A hypergraph allows to use edges that link more than just two compounds (see Figure 1.3b). A hypergraph H can be again denoted by a pair (V, E) where E is a set of hyperedges. While an edge is a set of two vertices, a hyperedge can be an arbitrary subset of vertices. However, this representation does not allow to distinguish between two sets of vertices that are necessary to represent the substrates and products of a reaction. Additionally, as mentioned before, metabolic networks can also contain reactions that are irreversible. To enable a more precise representation of the metabolic network, a directed hypergraph as described in [START_REF] Ausiello | Directed hypergraphs: Problems, algorithmic results, and a novel decremental approach[END_REF]) can be used. A directed hypergraph is a pair (V, A) where A is the set of directed hyperedges, that is the set of hyperarcs. An hyperarc e, is an ordered pair (t, h) of two sets of vertices that are called tail and head of e. They are denoted by tail(e) and head(e). Consequently, substrates and products of a reaction can be represented by these two sets. Here, tail(e) corresponds to the substrates and head(e) to the products of the reactions. Additionally, the distinct orientation of a hyperarc allows to represent the direction of the corresponding reaction (see Figure 1.3c). To [START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF]. (a) An undirected graph. The edges represent reactions that transform one compound into another one. It is difficult to represent reactions that have multiple substrates and/or products. This means that the three displayed edges A ↔ B, B ↔ C and B ↔ C could correspond to three different reactions but they could also be the result of, for example, the reaction B ↔ A + C + D. Since no distinct directions can be represented, it is not possible to distinguish between reversible and irreversible reactions. (b) An undirected hypergraph. A hyperedge has one arbitrary subset of vertices. Therefore, it is not possible to distinguish between the substrates and products of a reaction.

Here, the represented reaction could be B ↔ C + D but other interpretations like C ↔ B + D are also possible. Additionally, it is not possible to represent that a reaction is irreversible. (c) A directed hypergraph. The possible directions for each reaction are displayed by arrows. This representation allows for a simple discrimination between reversible and irreversible reactions. Additionally, there is a clear separation between substrates and products. Reaction A ↔ B is reversible whereas B → C +D is irreversible.

Stoichiometric matrix

All presented graph models represent the structure (topology) of the metabolic network but they do not take into account any information that might be available about the quantitative relationships between the substrates and the products of a reaction. Depending on the application, it might be important to include the stoichiometry of the reactions. A stoichiometric coefficient describes the quantity in which a metabolite participates in a specific reaction. For example, given the reaction 2A + 1B → 3C, two molecules of metabolite A and one molecule of metabolite B are needed to produce three molecules of C. It is possible to add the stoichiometric information as labels in a graph model. However, another commonly used representation is the stoichiometric matrix.

The stoichiometric matrix S is an m × n matrix where the m rows correspond to the metabolites in the metabolic network and the n columns to the reactions. The entry S ij corresponds to the stoichiometric coefficient of metabolite i in reaction j. If the stoichiometric coefficient is positive, the metabolite is produced by the corresponding reaction. If it is negative, the metabolite is consumed.

If the stoichiometric coefficient is zero, it indicates that the metabolite is not participating in this reaction or that the production compensates precisely its consumption.

The stoichiometric matrix itself does not contain any information about the reversibility of the reactions. Usually, this information is stored alongside the matrix in the form of lower bounds LB j and upper bounds U B j for the flux of each reaction j. For irreversible reactions, the lower bounds are greater or equal to zero. Changing the lower and upper bounds for source reactions can be used The network can also be represented as stoichiometric matrix S.

S 1 S 2 A C B T r 1 1 r 2 1 r 3 1 1 r 5 2 1 1 r 4 2 1 1 1 1 r 6 r 7 1 (a) S = r 1 r 2 r 3 r 4 r 5 r 6 r 7                     1 0 -1 0 0 0 0 S 1 0 1 0 0 -2 0 0 S 2 0 0 1 -2 -1 0 0 A 0 0 0 0 1 -1 0 C 0 0 0 0 0 1 0 B 0 0 0 1 0 1 -1 T (b)
The rows of S correspond to the metabolites in M and the columns to the reactions in R. The stoichiometric values are displayed as coefficients of the matrix. If a coefficient is positive, the corresponding metabolite is produced by the reactions. If it is negative, it is consumed. Hence, in the matrix representation, a direction is assumed for each reaction. However, it is not possible to see if a reaction is reversible or not. Lower and upper bounds for reactions have to be stored separately.

to model different growth media by allowing or restricting the uptake of specific metabolites, that is the import of metabolites from an external source.

Constrained-based modeling for metabolic networks

The stoichiometric matrix representation is used in constraint-based modeling of metabolic networks [START_REF] Covert | Constraints-based models: regulation of gene expression reduces the steady-state solution space[END_REF][START_REF] Palsson | The challenges of in silico biology[END_REF][START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF][START_REF] Lewis | Constraining the metabolic genotypephenotype relationship using a phylogeny of in silico methods[END_REF]. In constraintbased modeling, the idea is to analyze or identify possible flux distributions in the metabolic network under certain constraints.

The flux vector describes the state of all reactions of the metabolic network and shows if a reaction has a flux, i.e. is active. The flux vector v is an n-vector in which v j corresponds to the flux of reaction j.

A common assumption is that the metabolic network is in steady-state which means that every metabolite is produced in the same amount as it is consumed. Based on the stoichiometric matrix and the flux vector, these conditions can be formulated as constraints:

S • v = 0 (1.16) LB j ≤ v j ≤ U B j ∀j ∈ R rev (1.17) 0 ≤ v j ≤ U B j ∀j ∈ R irr .
(1.18) Equation 1.16 defines the mass balance for all metabolites. Additionally, lower and upper bounds are added to restrict the flux of single reactions (Equation 1.17). For irreversible reactions, the lower bound is set to zero (Equation 1.18). These constraints describe the core idea for several constraintbased approaches for metabolic networks. In the following, some of the most common among such approaches are presented that can also be used as a basis for more advanced methods. Probably the most important approach is Flux-balance-analysis (FBA) which is widely used to analyze metabolic networks (Orth et al., 2010b;[START_REF] Mahadevan | The effects of alternate optimal solutions in constraint-based genome-scale metabolic models[END_REF][START_REF] Edwards | Robustness analysis of the escherichiacoli metabolic network[END_REF][START_REF] Durot | Genome-scale models of bacterial metabolism: reconstruction and applications[END_REF][START_REF] Covert | Constraints-based models: regulation of gene expression reduces the steady-state solution space[END_REF][START_REF] Gottstein | Constraint-based stoichiometric modelling from single organisms to microbial communities[END_REF]. FBA uses the aforementioned constraints to restrict the feasible flux space. Additionally, FBA optimizes a specific objective function (Equation 1.19), e.g. it can simulate the growth of an organism by maximizing the production of biomass compounds.

max c T v (1.19) s.t. S • v = 0 (1.20) 0 ≤ v j ≤ U B j ∀j ∈ R irr (1.21) LB j ≤ v j ≤ U B j ∀j ∈ R rev (1.22)
Common applications for FBA are, for instance, to predict the growth of an organism on different media or under different conditions (e.g. aerobic and anaerobic). It can predict the phenotype [START_REF] Edwards | Robustness analysis of the escherichiacoli metabolic network[END_REF][START_REF] Edwards | In silico predictions of escherichia coli metabolic capabilities are consistent with experimental data[END_REF] and hence, it is also used to analyze the altered behavior of an organism that is subject to metabolic interventions, i.e. to predict the behavior of a mutant strain. The maximization of biomass or energy (ATP) production are frequently used objectives [START_REF] Feist | The biomass objective function[END_REF][START_REF] Pramanik | Stoichiometric model of escherichia coli metabolism: Incorporation of growth-rate dependent biomass composition and mechanistic energy requirements[END_REF]. The core idea of FBA is also used as basis for other similar methods, such as flux variability analysis. In Figure 1.5, a general work-flow of FBA is presented. A variant of FBA, parsimonious enzyme usage FBA (pFBA) [START_REF] Lewis | Omic data from evolved e. coli are consistent with computed optimal growth from genome-scale models[END_REF], assumes that an optimal solution should correspond the lowest overall flux through the network which simulates the minimization of the total amount of enzyme mass that is necessary. As argued in [START_REF] Segrè | Analysis of optimality in natural and perturbed metabolic networks[END_REF], when analysing metabolic networks, the maximization of the biomass production can be a reasonable assumption because it reflects the hypothesis that organism tend to optimize their growth under evolutionary pressure. Problematically, this assumption might not be valid for mutant strains that are the result of metabolic engineering because they are less exposed to evolutionary pressure than the original wild type strain.

Based on these observations, the authors in [START_REF] Segrè | Analysis of optimality in natural and perturbed metabolic networks[END_REF] present their method Minimization of metabolic adjustment (MOMA) which is similar to FBA but instead of maximizing the biomass to predict flux distributions for the mutant strain, they are minimizing the distance between the wild type flux and the flux in the mutant type (Equation 1.23). Here, the assumption is that the perturbed flux should remain as close as possible to the original flux because the organism might try to make as little changes as possible to its metabolic activity because using different reactions that were inactive before might also require different or additional enzymes.

min f (v) = |R| j=1 (v j -vj ) 2 (1.23)
MOMA is a quadratic optimization problem because the objective function minimizes the square root of the sum of the squared distances between the wild type flux v and the perturbed flux in the mutant type v. The flux values for the wild type can, for instance, be obtained by doing a FBA that maximizes the biomass production. In many cases, there is however not only one unique flux vector that leads to the optimal biomass production [START_REF] Mahadevan | The effects of alternate optimal solutions in constraint-based genome-scale metabolic models[END_REF]. This can render the choice for the wild type flux more difficult and the outcome of MOMA may differ depending on which exact flux vector was chosen as wild type. Alternative solutions are a general problem of FBA because it only computes one flux distribution. One possible approach to analyze alternative solutions is to identify by how much the flux for reactions can differ under the optimal condition, e.g. when the biomass is maximized. Flux variability analysis (FVA) allows to analyze the possible variation of the flux for all reactions under a specific condition [START_REF] Mahadevan | The effects of alternate optimal solutions in constraint-based genome-scale metabolic models[END_REF][START_REF] Burgard | Minimal reaction sets for escherichia coli metabolism under different growth requirements and uptake environments[END_REF]. It determines the minimum and the maximum flux values for all reactions under the optimal condition (Equation 1.24). This is achieved by introducing a constraint that fixes the biomass production to the optimal value z * calculated by FBA (Equation 1.28).

min/max v j (1.24)

s.t. S • v = 0 (1.25) 0 ≤ v j ≤ U B j ∀j ∈ R irr (1.26) LB j ≤ v j ≤ U B j ∀j ∈ R (1.27) c T v = z * (1.28)
In [START_REF] Lee | Recursive milp model for finding all the alternate optima in lp models for metabolic networks[END_REF], a recursive MILP is proposed to find alternative solutions. Furthermore, in [START_REF] Kelk | Optimal flux spaces of genomescale stoichiometric models are determined by a few subnetworks[END_REF] it is described that often, the possible flux distributions of a few subnetworks lead to the combinatorial explosion of optimal fluxes for the whole network.

A similar idea to FVA can also be used to analyze how reactions depend on each other if reactions are blocked under certain conditions [START_REF] Burgard | Flux coupling analysis of genome-scale metabolic network reconstructions[END_REF]. Flux coupling analysis identifies different relations between two fluxes v 1 and v 2 : They are directionally coupled (v 1 → v 2 ) if a non-zero flux for v 1 implies that v 2 must have a non-zero flux. This does not have to imply that a non-zero flux for v 2 always leads to a non-zero flux for v 1 . Furthermore, they are partially coupled (v 1 ↔ v 2 ) if a non-zero flux for v 1 implies a non-zero flux for v 2 and a non-zero flux for v 2 implies a non-zero flux for v 1 . Finally, two fluxes are fully coupled (v 1 ⇔ v 2 ) if the same conditions hold as for the partial coupling and in addition the implied fluxes are not variable but fixed which means that the ratio between the two fluxes is constant. If none of these relations holds for a pair of reactions, they are uncoupled. The differences are also displayed in Figure 1.6.

In [START_REF] Burgard | Flux coupling analysis of genome-scale metabolic network reconstructions[END_REF], the authors also identify blocked reactions which are reactions that can only have a zero flux under the given constraints that model, for instance, a certain medium, thereby restricting certain uptake fluxes (Equation 1.31). All reversible reactions are split into two irreversible reactions. Thus, negative fluxes are not possible (Equation 1.32) and it is sufficient to maximize the flux for a reaction to determine if the reaction is blocked or not. If reversible reactions are not split, both the maximum and minimum have to be computed which shows the direct relation to the FVA.

max v j (1.29) s.t. S • v = 0 (1.30) v uptake j ≤ v uptake max j ∀j ∈ R transport (1.31) v j ≥ 0 ∀j ∈ R (1.32)
Identifying blocked reactions is of interest because it can help to simplify the underlying network structure for specific cases. If reactions are blocked under the conditions that are modeled, it is possible to remove these reactions for other analyses under the same conditions. However, it is important to note that whether a reaction is blocked or not is very depending on the specific conditions (e.g. steady-state or uptake rates).

Another way to analyze the steady-state flux distributions in metabolic networks are elementary flux modes (EFM) [START_REF] Schuster | On elementary flux modes in biochemical reaction systems at steady state[END_REF]. A more recent review on the calculation and application of EFMs is given in [START_REF] Zanghellini | Elementary flux modes in a nutshell: properties, calculation and applications[END_REF]. The support of a flux vector v is defined by supp(v) = {i|v i = 0} that is the set of indices of non-zero elements. A flux vector v is called mode if it is non-trivial (i.e. v = 0) and is feasible under steady-state conditions. An EFM is a mode whose support cannot be written as proper subset of any other feasible mode. Therefore, EFMs are non-decomposable and represent minimal functional building blocks. A closely related idea are extreme pathways which are a subset of elementary modes [START_REF] Schilling | Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective[END_REF].

Predicting optimal intervention strategies

One of the areas of synthetic biology is the genetic engineering of microorganisms to improve the production of important target chemicals by a microorganism. Microorganisms are able to produce chemical compounds that can, for instance, be used in industry. The desired compound might often just be a by-product of the usual metabolism when the microorganism is growing and the yield can be very low. Thus, it is necessary to modify the metabolism of the organism to increase the yield of the desired target compound and to establish a more efficient production. The idea is to deactivate (knockout) specific reactions by manipulating the metabolic genes that are coding for the catalyzing enzymes. Depending on the role the knocked out reactions play in the metabolism, the organism might use different pathways to compensate for its lack.

Given the complexity of the problem, there is a need to identify the most promising intervention strategies in silico and it is necessary to predict and analyze the altered behavior after a knockout is introduced. A review on constraint-based methods for optimal strain design in silico is given in [START_REF] Maia | In silico constraint-based strain optimization methods: the quest for optimal cell factories[END_REF]. In the following, I will summarize some computational approaches to predict optimal intervention strategies.

In [START_REF] Burgard | OptKnock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization[END_REF], the authors present OptKnock which is one of the most well known methods for identifying knockouts. Their method is based on a bi-level optimization problem that aims to find a flux distribution that maximizes the production of the compound of interest while also maximizing the production of the biomass given some knocked out reactions.

max yj v target (1.33) s.t.        max vj v biomass s.t. S • v = 0 LB j (1 -y j ) ≤ v j ≤ U B j (1 -y j ) ∀j ∈ R        (1.34) j∈R y j = K (1.35) y j ∈ {0, 1}, v j ∈ R ∀j ∈ R (1.36)
The outer problem of the bi-level framework maximizes the bioengineering objective, i.e. the production of the target (Equation 1.33). The inner problem maximizes a cellular objective like the biomass production (Equation 1.34). The knockouts are modeled using a binary variable y j for each reaction 1 Theoretical background j. If the binary is set to one, the lower and upper bounds of the corresponding reaction are forced to zero which imitates a knockout (Equation 1.34). The approach identifies the optimal knockout set for a maximum number of knockouts (Equation 1.35). OptKnock is integrated in the COBRA Toolbox [START_REF] Heirendt | Creation and analysis of biochemical constraint-based models using the cobra toolbox v. 3.0[END_REF] available in MATLAB.

OptKnock can be overoptimistic because it assumes that both biomass and target production are maximized by the organism. For a given biomass production, even though a specific maximum target yield is possible, it does not mean that it will be how the organism behaves in reality. As explained in Section 1.2.4, there can be lots of alternative flux distributions that lead to the same optimal value. To account for competing pathways, in [START_REF] Tepper | Predicting metabolic engineering knockout strategies for chemical production: accounting for competing pathways[END_REF], the authors present another bi-level program which is called RobustKnock. It is similar to OptKnock. However, the crucial difference is that in the outer problem the minimum of the bioengineering objective is maximized. In this way, they ensure that the target metabolite has to be produced at least in the determined minimum amount in all flux distributions that are still achievable after knockouts have been introduced.

The approaches considered so far only model knockouts. To consider also up-an downregulation of fluxes, OptForce is presented in [START_REF] Ranganathan | OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions[END_REF]. Like OptKnock, it is a bi-level optimization problem. In a first step, for each reaction, the possible lower and upper bounds are computed by FVA to characterize the wild-type. This information is used to identify reactions sets that must change to achieve a desired overproduction specified by the user. These sets are called MUST sets. Based on the them, FORCE sets are computed which are minimal sets of reactions that must be genetically manipulated to force a change and to achieve the desired production yield.

Another different concept that can be used to identify intervention strategies in metabolic networks are minimal cut sets (MCS) which are described in [START_REF] Klamt | Minimal cut sets in biochemical reaction networks[END_REF][START_REF] Klamt | Generalized concept of minimal cut sets in biochemical networks[END_REF][START_REF] Ballerstein | Minimal cut sets in a metabolic network are elementary modes in a dual network[END_REF]. In general, MCSs are suitable to determine reactions that can be used to block metabolic functionalities. With regard to a specific target reaction that should be blocked, a set of reactions is a cut set if after all these reactions have been removed (e.g. knocked out), there is no feasible flux distribution possible that contains the target reaction. A cut set C is minimal if no proper subset of C is also a cut set.

The idea of MCS led to the introduction of constraint minimal cut sets (cMCS) in [START_REF] Hädicke | Computing complex metabolic intervention strategies using constrained minimal cut sets[END_REF]. Whereas MCSs are used to block undesired phenotypes, cMCSs allow to specify desired behaviors that should be kept. It prevents that MCSs are computed that block important functions. This concept was further generalized to regulatory MCSs that are able to not only model the knockout of reactions but also up/downregulations of reaction rates [START_REF] Mahadevan | Genome-scale strain designs based on regulatory minimal cut sets[END_REF]. Regulatory MCSs and cMCSs were already applied successfully to identify intervention strategies in genome-scale metabolic networks (von Kamp and Klamt, 2014;[START_REF] Mahadevan | Genome-scale strain designs based on regulatory minimal cut sets[END_REF].

Metabolic shifts

Metabolomics is a field that concerns itself with the measurement of metabolites [START_REF] Roessner | What is metabolomics all about?[END_REF]. It studies the metabolome which refers to the total number of metabolites that are found in the cells of an organism. Experiments to measure metabolite levels can be targeted, which means that a specific set of metabolites is quantified that are suspected to be relevant for the given experiment but there are also global analyses techniques available that take into account all measurable metabolites [START_REF] Fiehn | Metabolomics-the link between genotypes and phenotypes[END_REF].

Like transcriptomics and proteomics, metabolomics plays an important part in understanding how an organism reacts to changes in the environment. This can be, for example, the adaption of the organism to a change in the medium (e.g. nutritional stress). Furthermore, also the impact of genetic perturbations (e.g. mutations) on the organism can be analysed. Metabolomics makes it possible to explain how the metabolic activity of an organism changes after a perturbation and to identify which parts of the organism are affected. The change in metabolic activity is called metabolic shift.

Understanding metabolic changes under different conditions is important for advancing research in different fields like bioengineering or human health [START_REF] Sevin | Biological insights through nontargeted metabolomics[END_REF][START_REF] Goel | Metabolic shifts: a fitness perspective for microbial cell factories[END_REF].

In enrichment analyses of metabolomic data based on the metabolites and their relative abundances, pathways are identified that are more impacted by the perturbation than others [START_REF] Booth | Computational tools for the secondary analysis of metabolomics experiments[END_REF]. Several tools are already published that facilitate the analysis of metabolomic data. For example, MetaboAnalyst [START_REF] Xia | Metaboanalyst 3.0-making metabolomics more meaningful[END_REF][START_REF] Chong | MetaboAnalyst 4.0: towards more transparent and integrative metabolomics analysis[END_REF]) is a web-based tool that includes different modules for the integration of metabolomic data. Amongst others, it provides modules for metabolite set enrichment analysis (Xia and Wishart, 2010b), metabolic pathway analysis (Xia and Wishart, 2010a), as well as two-factor and time-series analyses [START_REF] Xia | Metatt: a web-based metabolomics tool for analyzing time-series and two-factor datasets[END_REF]. MetExplore [START_REF] Cottret | Metexplore: collaborative edition and exploration of metabolic networks[END_REF], another web-based application, is a versatile tool for the analysis of metabolic networks. It also implements MetaboRank [START_REF] Frainay | Metaborank: network-based recommendation system to interpret and enrich metabolomics results[END_REF] which can be used to interpret and enrich metabolomic data. MBRole [START_REF] Chagoyen | Mbrole: enrichment analysis of metabolomic data[END_REF], which is another web-server, can also be used to perform enrichment analysis of metabolomic data. Different tools and methods for the enrichment analysis of metabolomics are compared and evaluated in [START_REF] Marco-Ramell | Evaluation and comparison of bioinformatic tools for the enrichment analysis of metabolomics data[END_REF][START_REF] Booth | Computational tools for the secondary analysis of metabolomics experiments[END_REF][START_REF] Alonso | Analytical methods in untargeted metabolomics: state of the art in 2015[END_REF].

As also mentioned in [START_REF] Frainay | Computational methods to identify metabolic sub-networks based on metabolomic profiles[END_REF], one of the main issues in the analysis of metabolomic data is that there is no single metabolomics technology that allows to measure all metabolites. This means that only a partial view of the metabolome can be obtained which makes the subsequent analysis more complicated. Furthermore, as summarized in [START_REF] Booth | Computational tools for the secondary analysis of metabolomics experiments[END_REF] there is still a large number of metabolites that is unidentified and therefore also not characterized or annotated in known databases. A lack of refined networks and databases for organisms that are less commonly studied makes the analysis of metabolomic data more challenging.

Introduction

Microorganisms are already used on an industrial scale to produce important target chemicals. Metabolic engineering of strains has been widely used in order to optimize the bioconversion pathways aiming at higher product yield. However, accumulation of the target chemical can often negatively impact the cultivation of the microorganism [START_REF] Mukhopadhyay | Tolerance engineering in bacteria for the production of advanced biofuels and chemicals[END_REF]. Problematically, there are even more extreme cases where the target chemical is toxic for the used microorganism which might therefore restrict the efficiency of the production. Strain engineering also focuses on alleviating these bottlenecks, reducing cellular burdens that might limit product yield in order to achieve optimal production. Among such chemicals, there is a growing interest in ethanol due to its use as bio-fuel [START_REF] Mussatto | Technological trends, global market, and challenges of bio-ethanol production[END_REF]. Yeasts can convert sugars to ethanol during a process called fermentation [START_REF] Lin | Ethanol fermentation from biomass resources: current state and prospects[END_REF][START_REF] Bai | Ethanol fermentation technologies from sugar and starch feedstocks[END_REF][START_REF] Boulton | Yeast and biochemistry of ethanol fermentation[END_REF] which is already used to produce ethanol industrially.

The production of ethanol in yeast can be improved by metabolic engineering [START_REF] Borodina | Advances in metabolic engineering of yeast saccharomyces cerevisiae for production of chemicals[END_REF][START_REF] Nielsen | Metabolic engineering of yeast for production of fuels and chemicals[END_REF][START_REF] Van Vleet | Yeast metabolic engineering for hemicellulosic ethanol production[END_REF]. However, an important limiting factor for the production of ethanol in yeast on an industrial scale is indeed the toxicity of ethanol for the yeast.

When ethanol accumulates, it inhibits growth but leads also to a decline of the ethanol production [START_REF] Van Uden | Ethanol toxicity and ethanol tolerance in yeasts[END_REF][START_REF] Dombek | Magnesium limitation and its role in apparent toxicity of ethanol during yeast fermentation[END_REF][START_REF] D'amore | Ethanol tolerance of yeast[END_REF][START_REF] Casey | Ethanol tolerance in yeasts[END_REF][START_REF] Lam | Engineering alcohol tolerance in yeast[END_REF][START_REF] Panchal | A study of ethanol tolerance in yeast[END_REF]. Therefore, different approaches are investigated to increase the tolerance of yeast to ethanol [START_REF] Dombek | Magnesium limitation and its role in apparent toxicity of ethanol during yeast fermentation[END_REF][START_REF] Alper | Engineering yeast transcription machinery for improved ethanol tolerance and production[END_REF][START_REF] Lam | Engineering alcohol tolerance in yeast[END_REF]. Such tolerance can differ between yeast strains [START_REF] Casey | Ethanol tolerance in yeasts[END_REF].

As mentioned in Section 1.3, different constraint-based approaches were already developed that can be used to identify optimal intervention strategies in metabolic networks. The methods that were presented (e.g. OptKnock, RobustKnock) are some of the most common ones for metabolic engineering but they are mostly based on single-objective linear programs or bi-level approaches. However, in general, multi-objective optimization is already commonly applied in bioinformatics and computational biology [START_REF] Handl | Multiobjective optimization in bioinformatics and computational biology[END_REF]. Multi-objective approaches have the advantage that they can investigate conflicting objectives at the same time.

In the context of metabolic engineering, this can be, for example, the maximization of a cellular objective (e.g. biomass production) and the maximization of an engineering objective (e.g. the production of a target chemical). For example, in [START_REF] Patané | Multiobjective optimization of genome-scale metabolic models: the case of ethanol production[END_REF], the authors propose a multi-objective metabolic engineering algorithm that can model gene knockouts but also up-and downregulation of enzymes. Their method solves a multi-objective optimization problem for biomass production and ethanol production. Other examples of the application of multi-objective frameworks for ethanol production are presented in [START_REF] Vera | Multicriteria optimization of biochemical systems by linear programming: application to production of ethanol by saccharomyces cerevisiae[END_REF][START_REF] Sendín | Model based optimization of biochemical systems using multiple objectives: a comparison of several solution strategies[END_REF]. In [START_REF] Andrade | Momo-multi-objective metabolic mixed integer optimization: application to yeast strain engineering[END_REF], a multi-objective approach called MOMO is introduced that identifies points in the Pareto front that represent different tradeoffs between biomass and target production. Afterwards, MOMO enumerates and analyzes possible single knockouts for each point in the Pareto front. These methods however do not account for the fact that ethanol is toxic for yeast.

In this chapter, we present a new constraint-based approach inspired by MOMO [START_REF] Andrade | Momo-multi-objective metabolic mixed integer optimization: application to yeast strain engineering[END_REF] that proposes knockout sets to improve the production of any target chemical that is also toxic for the microorganism and whose accumulation might inhibit growth and slow down the production. We therefore also ensure that the introduced knockouts do not restrict the metabolic activity of the organism that is crucial for a better resistance against the toxic target chemical. The approach is separated in two parts. In the first part, a multi-objective optimization problem is used to identify different efficient flux distributions that maximize the biomass production, the target production and the resistance against the toxic target. In the second part, different approaches are presented to enforce a flux distribution of interest that was identified in the first part. For the second part, different ideas were explored. As a first idea, a MILP is proposed to enumerate different knockouts. Due to its limitations, there was a need to develop a second idea which is based on identifying and cutting off subnetworks. The present method is applied to the case-study of ethanol production in yeast.

Computation of tradeoffs

In the first part of this approach, a multi-objective optimization problem is used to identify tradeoffs between biomass production, target production and the resistance against the toxic target. To calculate a score that describes the resistance of the organism against the toxic target, it is necessary to have information on certain critical reactions that must have been identified beforehand. A critical reaction is a reaction that increases the resistance of the organism against the toxic target when it is active. Computing these tradeoffs allows to gain an overview of optimal outcomes for which feasible flux distributions exist.

Toxicity resistance score

To evaluate the resulting resistance against the toxic target for a specific flux distribution, a score r is defined that depends on the flux values of the critical reactions:

r = j∈Rcritical (ω j t j + ω j v j U B j ). (2.1)
The score is based on three different assumptions concerning the reactions that are critical for the resistance against the target metabolite:

1. A non-zero flux is more important than the amount of flux for a critical reaction. This assumption is represented by the first part of the sum. For each reaction j that was identified to be important for the toxicity resistance, a binary variable t j is introduced. If the corresponding reaction has a flux, the binary variable is set to one. Thus, if many important reactions have a flux, the score will be high. The exact quantity of the flux has no impact on this part of the score.

2. A higher flux is better than a lower one. It is represented by the second part of the sum. The flux values v j are normalized by their upper bound U B j to be between 0 and 1 before they are added to the score to ensure that they are less significant for the total score than the binary variables.

3. Some reactions have a higher impact on the toxicity resistance than other reactions. These differences can be modeled by setting individual weights w j for each critical reaction.

Multi-objective optimization problem

The toxicity resistance score is used to formulate a multi-objective optimization problem that maximizes three different objectives: The biomass production v biomass , the target production v target and the toxicity resistance score r.

max v target , v biomass , j∈Rcritical (ω j t j + ω j v j U B j ) (2.2) s.t. S • v = 0 ∀i ∈ M (2.3) LB j ≤ v j ≤ U B j ∀j ∈ R (2.4) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.5) v j , ω j ∈ R; t j ∈ {0, 1} (2.6)
The system assumes a steady-state (Equation 2.3) and for each reaction upper and lower bounds are specified (Equation 2.4). Furthermore, constraints are introduced to connect the flux values v j to the binary variables t j for all critical reactions R critical (Equation 2.5). If the corresponding flux value is zero, the binary variable must be set to zero too. If the reaction has a non-zero flux, the binary variable can be set to one. Here, it is not necessary to explicitly model v j = 0 =⇒ t j = 1 that would ensure that a binary variable has to be set to one if the corresponding flux value is different from zero because the optimization problem is maximizing the toxicity resistance score. This means that if a t j can be set to one, the solver will set it to one to increase the score.

In practice, a small threshold m is used to identify reaction fluxes different from zero and the implication constraint is remodeled:

-∞ ≤ t j -v j ≤ 1 -m. (2.7)
There are two possibilities. If the reaction has a flux greater or equal to m, the difference d j = v j -m is greater to equal to zero. Therefore, t j ≤ 1 + d j holds and t j can be set to one. If the reaction flux is smaller than m (and hence it is assumed to be zero), d j is negative and t j can only be set to a value smaller than one, and since t j is a binary, it can only be set to zero.

However, this is only viable if the reaction cannot have a negative flux. For reversible reactions, the constraint needs to be adjusted:

-∞ ≤ t j -|v j |≤ 1 -m. (2.8)
To model the absolute value |v j |, this part can be split:

-∞ ≤ t j -v + j -v - j ≤ 1 -m.
(2.9)

In this case, v + and v -are two new non negative variables with the corresponding bounds:

0 ≤ v + j ≤ U B j • y j (2.10) 0 ≤ v - j ≤ |LB j |•(1 -y j ). (2.11)
The idea is the same as splitting a reversible reaction in two irreversible reactions. The forward direction is represented by v + j and the backward direction by v - j . Since in the case of a reversible reaction, the lower bound is negative, the absolute value of it is used as upper bound for v - j . The new binary variable y j prevents that v + j and v - j are non zero at the same time. The variables v + j and v - j are only needed if the reaction is reversible (i.e. the lower bound is negative). Otherwise, Equation 2.7 is sufficient.

In the resulting Pareto front of this multi-objective optimization problem, each point describes a different possible tradeoff between the three objectives. Computing only the extreme points is enough to gain a broad overview of the whole Pareto front. The idea is to choose a point in the Pareto front that has a high target yield and a sufficient biomass production and to find a way to enforce the values for target and biomass yield and the toxicity score. One possibility is to introduce knockouts that restrict the phenotypic space around the desired values for target and biomass production as mentioned in Section 1.3.

The multi-objective optimization problem was implemented in C++ and PolySCIP [START_REF] Borndörfer | Polyscip[END_REF] was used as solver which is part of the SCIP Optimization Suite 5.0 (Gleixner et al., 2017a,b). The metabolic network was modeled using the library metnetlib which is available on https://gitlab.inria.fr/erable/kirikomics/metnetlib.

2.3 First approach -Identifying knockouts that enforce the tradeoff flux

MILP formulation

After a promising point p of the Pareto front has been picked, a modified version of the first MILP is used to search for knockouts that might restrict the phenotypic space enough to enforce the desired values for the toxicity score and the target and biomass yield. The values are fixed to the corresponding values from the Pareto front v p target , v p biomass and r p (Equations 2.15, 2.16, 2.17). It is not necessary to minimize or maximize an objective function because we are looking for different feasible solutions that correspond to these fixed values. Additionally, new binary variables x j are introduced to model reaction knockouts. It is possible that not all reactions are potential candidates for a knockout. If it is known before that it is not viable to knockout a certain reaction in vivo, it is not included in R knockout . Each reaction j in R knockout has an associated binary variable x j . If the binary is set to one, the reaction flux is forced to zero (Equation 2.18). Furthermore, the number of total knockouts that are introduced in the metabolic network are fixed to a specified number K (Equation 2.19).

s.t. j∈R s ij v j = 0 ∀i ∈ M (2.12) LB j ≤ v j ≤ U B j ∀j ∈ R (2.13) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.14) v target = v p target (2.15) v biomass = v p biomass (2.16) r = r p
(2.17)

x j = 1 =⇒ v j = 0 ∀j ∈ R knockout (2.18) x j = K ∀j ∈ R knockout (2.19) v j , ω j ∈ R; t j , x j ∈ {0, 1}; K ∈ N + (2.20)
To model x j = 1 =⇒ v j = 0, in CPLEX, indicator constraints can be used. They are also called IfThen constraints. In SCIP, the included indicator constraint is modeled as follows:

z = 1 =⇒ ax ≤ 0. (2.21)
Here, z is a binary variable. To model the equality in Scip if the flux can be negative, the constraint in Equation 2.18 has to be split into two new indicator constraints:

x j = 1 =⇒ v j ≤ 0 (2.22) x j = 1 =⇒ -v j ≤ 0. (2.23)
The proposed MILP can be used to propose reaction knockouts that still allow the desired values for the three objectives. To enumerate different combination of knockouts, it has to be solved multiple times and after each iteration, the computed knockout set has to be excluded as solution by adding a corresponding constraint to the optimization problem. A knockout set P contains all variables x j that were set to one by the solver. By adding Equation 2.24 as constraint, it can be prevented that the exact same combination of x j is chosen again:

xj∈P x j ≤ K -1. (2.24)
The problem can be solved repeatedly until it becomes infeasible to enumerate all knockout sets that are possible for the fixed values for biomass and target production and the toxicity score. For a small K, i.e. K = 1 or K = 2, it is possible to enumerate all knockout sets even for larger networks. Depending on the size of the network and on the size of R knockout , for larger K, this task becomes more difficult. Therefore, reducing the size of R knockout can help simplify the problem.

Reducing the number of knockout candidates

Several steps can be taken to reduce the size of R knockout . First of all, a pre-selection should be done. For example, it can be difficult to knockout some transport reactions, given that a substrate can sometimes be transported by more than one system and also because transporters can be nonspecific. Exchange reactions are artificial reactions added to the networks to model boundaries of the organism. Furthermore, all genome-wide reconstructed models have reactions that do not have any associated genes. These reactions should be removed from R knockout because their knockout is not applicable in vivo.

As a next step, a certain aspect of Flux coupling analysis is used to determine groups of reactions that are knocked out simultaneously. We are interested in identifying such groups of reactions. For each two reactions i, j in a group, it holds that:

v j == 0 ⇐⇒ v i == 0. (2.25)
To achieve this, single knockouts for all reactions currently in R knockout are done and afterwards, FVA is used to analyze the possible flux values for the other reactions in R knockout . If a reaction can only have a zero flux, the other direction of the equivalence of Equation 2.25 is verified. As a result, reactions will be grouped together and only one representative for each group is present in R knockout . All others are removed. If the representative is chosen as knockout, afterwards, it is possible to replace it by other reactions from the same group without changing the result. This step also eliminates the need to simplify linear pathways. For instance, if there is a linear pathway like A → B → C → D, it is possible to reduce it to A → D thereby lowering the number of metabolites and reactions without changing the overall behavior of the metabolic network. By determining the coupled groups of reactions and only using one representative that can be knocked out, it is also possible to remove parts of linear pathways from R knockout without changing the network.

To further reduce the size of R knockout , it might be important to analyze the influence of single knockouts. All remaining reactions in R knockout are knocked out separately (Equation 2.29) and afterwards, FBA is used to maximize the biomass production. If it is zero, the candidate k can be removed from R knockout . Otherwise, the target production is maximized subsequently. Likewise, if the maximum target production is zero, the candidate k is removed.

max v biomass (or v target ) (2.26) s.t. j∈R s ij v j = 0 ∀i ∈ M (2.27) LB j ≤ v j ≤ U B j ∀j ∈ R (2.28) v k = 0 (2.29) v j ∈ R (2.30)
The reactions that are removed in this way from R knockout cannot be chosen by the MILP presented in Section 2.3.1 because they do not allow for a biomass or target production. However, removing them beforehand leads to less binary variables in the MILP because only for reactions in R knockout , an associated binary variable x j is introduced.

All the steps to reduce the size of R knockout that are presented until now are independent of the point in the Pareto front that was chosen. The next steps however will depend on the exact point in the Pareto front and have to be repeated if different points are investigated.

For a specific tradeoff p, it is possible to verify for each candidate k if there exists a flux distribution where reaction k can have a zero flux (i.e. it can be knocked out). The corresponding MILPs are shown below (Equation 2.31 to 2.38).

min/max v k (2.31) s.t. j∈R s ij v j = 0 ∀i ∈ M (2.32) LB j ≤ v j ≤ U B j ∀j ∈ R (2.33) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.34) v target = v p target (2.35) v biomass = v p biomass (2.36) r = r p (2.37) v j ∈ R; t j ∈ {0, 1} (2.38)
To model the specific tradeoff, the values for the biomass and target productions and the toxicity score are fixed to the values of the tradeoff. For irreversible reactions, it is sufficient to minimize v k . If the minimum is greater than zero, it has to be active in all flux distributions that are possible for the fixed values of the tradeoff and it cannot be knocked out. Thus, it can be removed from R knockout . For reversible reactions, it needs to be verified that the minimum is smaller or equal to zero and that the maximum is greater or equal to zero. Otherwise, it will be removed from R knockout .

Again, the removed candidates cannot be chosen by the MILP in Section 2.3.1 but their early removal will reduce the number of binary variables in the MILP.

In the next step, the idea is to compare the FVA values between the wild type flux and the tradeoff flux for all remaining candidates in R knockout . Therefore, first, FVA is done for all reactions in R knockout for the wild type. Afterwards, the values for biomass production, target production and the toxicity resistance score are fixed to the values of the tradeoff and the FVA is repeated. Based on the comparison of the possible flux values between wild type and tradeoff flux, the reactions are sorted into groups. The idea is that some reactions might be more likely to have an impact on the resulting flux than others. It might be advantageous to try candidates that are more likely to change the flux first or, when looking for double (or larger) knockouts, to ensure that every knockout set contains at least one of the higher priority candidates.

The first group contains all reactions that must have a non zero flux in the wild type and that must have a zero flux in the tradeoff flux. The flux of these reactions must change to establish the tradeoff flux and therefore, these are potentially very strong candidates and they should have a high priority.

The assumption for the next group is similar but less strict. It includes all reactions that must have a non zero flux in the wild type and that can have a zero flux in the tradeoff. All remaining candidates are put into the last group.

Evaluation of knockout sets

After reducing the size of R knockout as much as possible, the MILP in Section 2.3.1 is used to compute different knockout sets that are possible for the chosen point in the Pareto front. The disadvantage of the given formulation is that the proposed knockout sets do not necessarily enforce the desired values for biomass production, target production and the toxicity resistance score but all knockout sets are computed that still allow for these values but that can also result in less optimal ones. Hence, it is necessary to evaluate all computed knockout sets afterwards to identify the ones that lead to the best results.

To evaluate a knockout set P , first, the biomass production is maximized after all reactions in P are knocked out.

max v biomass (2.39) s.t. j∈R s ij v j = 0 ∀i ∈ M (2.40) LB j ≤ v j ≤ U B j ∀j ∈ R (2.41) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.42) x j = 1 =⇒ v j = 0 ∀j ∈ R knockout (2.43) x j = 1 ∀j ∈ P (2.44) v j ∈ R; t j , x j ∈ {0, 1} (
v biomass = v * biomass (2.48) v j ∈ R; t j , x j ∈ {0, 1} (2.49)
Here, the target production is minimized to account for the worst case in which the organism does not prioritize the target production at all. As described in Section 1.3, due to alternative pathways, it is not enough to maximize the target production because this might not be in accordance with the behavior of the organism and no real improvement for the target yield might be achieved which also happened in some cases for OptKnock. The evaluation is based on the assumption that the biomass production is still maximized after introducing knockouts. As described in Section 1.2.4, this is not necessarily the case and it is also possible to use an approach like MOMA that minimizes the distance to the wild type flux. The toxicity resistance score however is maximized because if the organism is under the stress of a toxic environment, it should have a biological incentive to increase its own resistances.

Main drawbacks and other approaches

In the results that are presented in Section 2.5.3, it will become clear that this approach is not viable for larger networks. The required number of knockouts to enforce the desired tradeoff values might be too large and enumerating all possible combinations would take too much time and is not reasonable. After evaluating the computed knockout sets, it becomes clear that almost all of the enumerated knockout sets do not lead to a change in the target production compared to the wild type. The MILP is not restrictive enough. It is possible to enumerate single or double knockouts but already knockout sets of size three are problematic because there are too many combinations that are possible. Hence, it was necessary to change the approach.

A first idea was to address the main problem of the used MILP and render it more restrictive. This could be done by using a max-min problem (similar to RobustKnock) while computing the knockout sets. The problem should maximize the minimum ethanol production. In this way, it should be ensured that the introduced knockouts lead to a change in the target production.

Although the idea itself is simple, the main drawback is that the modified problem is no longer a linear optimization problem and solving it becomes considerably more difficult.

2.4 Second approach -Isolating the active subnetwork Metabolites that are exclusively connected to active reactions are referred to as internal metabolites.

Metabolites that are connected to at least one active reaction and at least one reaction with a zero flux are metabolites lying on the border of the hyperpath.

Furthermore, inactive reactions that are connected to at least one metabolite lying on the border of the hyperpath are called incoming or outgoing reactions. An inactive reaction is an incoming (outgoing) reaction if it is producing (consuming) a metabolite on the border. It is possible for a reaction to be incoming and outgoing at the same time.

Enumerating different hyperpaths

To enumerate topologically different solutions for one point p of the Pareto front, a MILP is solved. Again, we assume that the system is in steady-state (Equation 2.50), each reaction has lower and upper bounds (Equation 2.51) and the critical reactions are associated with binary variables (Equation 2.52). Additionally, the values for target production, biomass production and the toxicity score are fixed at their optimal values in p (Equation 2.53, 2.54, 2.55). Since we are interested in topologically different solutions, each reaction j is associated with a binary variable a j that indicates if the corresponding reaction has a non-zero flux and therefore participates in the solution (Equation 2.56). Additionally, smaller hyperpaths (solutions that have less active reactions) are preferable because fewer reactions of the whole network are needed. However, minimizing the number of active reactions means that j∈R a j has to be minimized which renders the problem computationally more expensive. Hence, only the size of the resulting hyperpath is limited to K (Equation 2.57) which does not compute the smallest hyperpaths but avoids that large hyperpaths appear as a solution.

s.t. j∈R s ij v j = 0 ∀i ∈ M (2.50) LB j ≤ v j ≤ U B j ∀j ∈ R (2.51) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.52) v target = v p target (2.53) v biomass = v p biomass (2.54) r = r p (2.55) a j = 0 ⇐⇒ v j = 0 ∀j ∈ R (2.56) j∈R a j = K (2.57) v j , ω j ∈ R; t j , a j ∈ {0, 1} (2.58)
After one hyperpath H is computed, it has to be excluded as solution by adding another constraint before the problem is solved again:

j∈H a j ≤ |H|. (2.59)
This process can be repeated until the problem becomes infeasible and therefore no more new topologically different solutions exist or until the desired number of hyperpaths have been enumerated.

Rating a hyperpath

After having enumerated different hyperpaths for one point in the Pareto front, it must be investigated which of these hyperpaths are preferable compared to the others. The idea is to compare what kind of production values a specific hyerpath can achieve. One hyperpath is representing a certain subnetwork which must be cut off from the rest of the metabolic network.

The hyperpath can be isolated by knocking out all incoming and outgoing reactions. Alternatively, it is possible to remove all sources that are not part of the hyperpath and knock out only the outgoing reactions. If all sources of the network are part of the hyperpath, it is enough to knock out all outgoing reactions. Indeed, doing this prevents that the flux can deviate from the hyperpath. Removing all external sources is equal to knocking out all incoming reactions because if all outgoing reactions are removed, incoming reactions can only have a flux if they are fed from an external source that is not part of the subnetwork. An illustration is shown in Figure 2.2.

With the hyperpath isolated from the rest of the network, it is now possible to compute the biomass production, target production and the toxicity score just for the remaining subnetwork. We assume that the organism is still maximizing the biomass production after knocking out all outgoing reactions. This assumption is not accurate in all cases but it is possible to change it if a better estimation is on hand and to modify the MILP accordingly.

The following MILP can be used to compute the maximum biomass production v * biomass in the hyperpath.

max v biomass (2.60) s.t. j∈R s ij v j = 0 ∀i ∈ M (2.61) LB j ≤ v j ≤ U B j ∀j ∈ R (2.62) v j = 0 =⇒ t j = 0 ∀j ∈ R critical (2.63) v j = 0 ∀j ∈ R incoming (2.64) v j = 0 ∀j ∈ R outgoing (2.65) v j , ω j ∈ R; t j ∈ {0, 1} (2.66)
As already mentioned, it is possible that R incoming = ∅.

Afterwards, minimum and maximum target productions and minimum and maximum toxicity resistance scores are calculated when the biomass production is fixed to its optimum v * biomass .

min/max v target (2.67)

s.t. (2.61) -(2.65) (2.68) v biomass = v * biomass (2.69) v j , ω j ∈ R; t j ∈ {0, 1} (2.70)
T B A C S 1 S 2 r 1 r 3 r 2 r 4 r 5 Figure 2.2:
Illustration how incoming reactions can be blocked by removing external sources. In this small example, the hyperpath is highlighted by the rectangle which means that r 1 and r 4 are the active reactions of the hyperpath. r 5 is an incoming reaction and r 2 is an outgoing reactions. The metabolites S 1 and S 2 are sources of the network. T is an internal metabolite. A, B and S 1 are metabolites on the border of the hyperpath. We are interested in producing the target metabolite T while also maintaining some biomass production (metabolite B). To make sure that the flux cannot deviate from the desired hyperpath flux distribution, all incoming and outgoing reactions have to be knocked out. An alternative way is to only knock out all outgoing reactions and to remove all external sources that remain outside of the hyperpath. After knocking out all outgoing reactions, the only way that the substrates of incoming reactions are present is through external sources of the network. So in this case, by knocking out reaction r 2 and removing the metabolite S 2 from the medium, C cannot be produced and the incoming reaction r 5 cannot take place.

v biomass = v * biomass (2.73) v j , ω j ∈ R; t j ∈ {0, 1} (2.74)
The result can be used to compare different hyperpaths. A hyperpath is rated high if it has a good target production. A high maximum target production is not sufficient. To account for the worst case, it is more important to have a high minimum target production. Preferably, minimum and maximum target productions are very close to each other. Furthermore, it should be possible to achieve high values for the toxicity resistance score. Contrary to the target production, for the toxicity resistance score, it is more important to have a high maximum value. Since we assume that the critical reactions might only be active if the need arises, that is when the concentration of the toxic product is high, a minimum score close to zero should not be a criterion to dismiss the corresponding hyperpath. In this case, it is more important to be able to achieve a high maximum toxicity resistance score because then, it is possible that a significant number of critical reactions can be activated if the organism is exposed to the toxic product.

Based on these assumptions, hyperpaths with high minimum target production and high maximum resistance score should be considered for further investigation. Hyperpaths that do not fit these criteria are less efficient.

Identifying smaller knockout sets

For simplification, the next steps assume that all sources of the network are included in the hyperpath. If this is not the case, it is necessary to remove the remaining sources (e.g. by removing them from the media). If there are sources outside of the hyperpath and it is not possible to remove these sources, not only the outgoing but also the incoming reactions have to be considered in the subsequent approaches. So far, to isolate the hyperpath from the rest of the network, all outgoing reactions are knocked out. However, the number of such reactions can be high and it might not be feasible to knockout such a high amount of reactions in vivo. Consequently, it is necessary to reduce the number of reactions that have to be knocked out. Our first ideas were based on minimal cut sets and topological precursors.

In the final approach, we decided to use a random exploration approach to find smaller subsets that are still leading to good values for minimum target production and maximum resistance score.

Topological precursor cut sets

The first approach is based on the identification of minimal topological sources for the incoming and outgoing reactions. The concept for topological precursors is introduced in [START_REF] Cottret | Enumerating precursor sets of target metabolites in a metabolic network[END_REF]Acuña et al., 2012b).

In (Acuña et al., 2012b), the notion of forward propagation is used to define what a precursor is. If M is a set of metabolites, then the forward propagation of M , denoted by F wd(M ), is the set of metabolites that can be produced from M . Subs(r) are all metabolites that are substrates of reaction r. Likewise, P rod(r) contains all metabolites that are products of reaction r. Reac(M ) is the set of reactions that can take place if all the metabolites in M are present in the network because

}) = {S 1 , S 2 , S 3 , A, C, T }.
all substrates of these reactions are available which means that Reac(M ) = {r ∈ R|Subs(r) ⊆ M }. Furthermore, if R is a set of reactions, the authors define two more sets: Subs(R) = ∪ r∈R Subs(r) and P rod(R) = ∪ r∈R P rod(r). Based on these notations, F wd(M ) can be computed by the recursion M i+1 = M ∪ P rod(Reac(M i )). The recursion starts from M 0 = M and finishes when a fixed point is reached which means that no more new metabolites can be added. An example is shown in Figure 2.3. Moreover, as proposed in [START_REF] Cottret | Enumerating precursor sets of target metabolites in a metabolic network[END_REF], the authors include other metabolites in their model which are called internal supply and are always available. F wd Z (M ) is the forward propagation of M given that the set of metabolites Z is an internal supply and the recursion can be reformulated as

M i+1 = M ∪ P rod(Reac(M i ∪ Z)).
The set of metabolites S denotes the set of all source metabolites of the metabolic network which means that they can be available as an external supply. Based on the idea of forward propagation and internal supply, a set of source metabolites X ⊆ S is a precursor set of the set of target metabolites T if F wd Z (X) contains both T and also Z. It is necessary that Z is included to make sure that it can be reproduced.

Since we are not interested in producing a set of target metabolites but our goal is to cut off all incoming and outgoing reactions, the most important idea is the precursor cut sets. The authors of [START_REF] Cottret | Enumerating precursor sets of target metabolites in a metabolic network[END_REF] define that a set of sources X ⊆ S is a cut set of the set of target metabolites T if and only if the set S \ X is not a precursor set of T . Hence, the production of T is prevented by removing the the elements in X as source metabolites. To ensure that the production of all metabolites in T is cut off, it is necessary to introduce a special target metabolite t * . For each metabolite t i ∈ T , a new reaction r i is added with Subs(r i ) = {t i } and P rod(r i ) = {t * }. Afterwards, the new set of target metabolites is T = {t * }.

Since we are interested in cutting off reactions, we have to modify the network because the idea of precursor cut sets is used to prevent the production of metabolites. One possibility is to add a new metabolite for each incoming and outgoing reaction that splits the reactions in two. For each incoming or outgoing reaction r i , a new metabolite m i is created and r i is split into r i and r i whereas Subs(r i ) = Subs(r i ), P rod(r i ) = {m i } and Subs(r i ) = {m i }, P rod(r i ) = P rod(r i ).

Furthermore, as explained before, a special target metabolite t * has to be created and one additional reaction for each m i must be added that consumes m i and produces t * .

To identify precursor cut sets for the outgoing reactions, the directions for all reactions of the network have to be reversed. This means that all metabolites that were sources of the network become sinks of the network and likewise, all metabolites that were previously sinks are now sources of the metabolic network. However, this implies that the precursor cut sets for incoming and outgoing reactions have to be identified separately. Furthermore, only metabolites that are not part of the hyperpath can be considered as sources because otherwise they are needed to obtain the biomass and target production. It might be preferable to cut out the identified hyperpath and to only take into account the remaining part of the network when looking for precursor cut sets.

As presented in (Acuña et al., 2012b), a minimal precursor cut set can be found quite easily. Starting from X = ∅, sources are added if the target cannot be produced. As a result, the set X = S \ X is a minimal precursor cut set. Identifying minimum precursor cut sets is however considerably more difficult but also more important because we are interested in the smallest changes that are necessary to achieve our goal. The authors in [START_REF] Cottret | Enumerating precursor sets of target metabolites in a metabolic network[END_REF] present an approach to enumerate all precursor sets that uses a kind of depth-first search on the hyperpath that traverses the reactions in their opposite direction to explore paths from the target to the sources. To avoid fake solutions, some cycles have to be removed from the metabolic network in a preprocessing step because they can lead to the production of certain metabolites without any external sources. During the preprocessing, these undesired cycles are removed by deleting some reactions. It is a heuristic and the preprocessed network is not always the same which also influences the minimal precursor sets that are computed afterwards.

In practice, we abandoned the idea to use topological precursors cut sets as a solution to our problem due to difficulties that were linked to the removal of the aforementioned cycles. Additionally, it can only identify sources as cut points and we are also interested in proposing reaction knockouts as intervention strategy. Furthermore, since this approach is based on the topology of the network, it might have been preferable to also remove small metabolites and co-factors from the network to simplify the network structure. It can however be difficult to automatically identify and remove cofactors in a way that the stoichiometry of the network is kept intact. The stoichiometry is not needed to identify minimal precursor sets but it might be important for subsequent analyses since in general, our method should include the stoichiometry of the network. It might be possible to circumvent some of the difficulties by using the idea of stoichiometric precursor sets [START_REF] Andrade | Enumeration of minimal stoichiometric precursor sets in metabolic networks[END_REF] but so far, we did not explore this concept further because we assumed that the computation might be more difficult and there might be other approaches that are more accessible for our problem.

Minimal cut sets

Next, we were interested in applying the concept of minimal cut sets to our problem. In this case, we want to find the smallest MCSs that block all outgoing reactions for a hyperpath. Therefore, all outgoing reactions are set as the target reactions. Additionally, it is not possible to knockout the active reactions of the hyperpath. The MSCEnumerator presented in (von Kamp and Klamt, 2014) enumerates MCSs and cMCSs.

An implementation is available in CellNetAnalyzer [START_REF] Klamt | An application programming interface for cellnetanalyzer[END_REF][START_REF] Klamt | Structural and functional analysis of cellular networks with cellnetanalyzer[END_REF][START_REF] Von Kamp | Use of cellnetanalyzer in biotechnology and metabolic engineering[END_REF]. CellNetAnalyzer uses MATLAB.

In (von Kamp and Klamt, 2014), MSCEnumerator is used to compute the smallest MCSs in a genome-scale metabolic network. However, as shown in their results, it can be computationally expensive for larger networks to compute MCSs. In our case, this approach was not suitable in practice. We assume that one of the main problems is that a larger number of reactions has to be knocked out to block all of the outgoing reactions.

Random exploration

In this last approach, we propose a random exploration that identifies smaller subsets of the outgoing reactions that are still leading to good values for minimum target production and maximum resistance score. The starting point for this approach are all outgoing reactions for a hyperpath that was chosen based on its promising values. The outgoing reactions are forming the reaction pool P 0 . The maximum biomass for this hyperpath is used as reference value for the biomass production v ref biomass .

Step 1: The reaction pool is split into two different subsets. For each reaction r in the reaction pool P 0 , all reactions in P 0 are knocked out except reaction r. Afterwards, the maximum biomass is computed. If the maximum biomass remains the same as v ref biomass , not knocking out only reaction r does not change the result and reaction r will be part of subset 2. If the maximum biomass is not the same as v ref biomass , not knocking out reaction r does lead to a change and reaction r will be part of the subset 1. This procedure is repeated for all reactions in P 0 . Consequently, subset 1 contains all reactions that have to be knocked out at all time because otherwise it is not possible to obtain the same production values as in 2.4.1. However, subset 2 contains reactions that do not have to be necessarily knocked out and it might be possible to remove several of them from the knockout set without changing the production values.

Step 2: A new knockout P 1 set is generated. All reactions of subset 1 must be in P 1 . One by one, reactions from subset 2 are randomly drawn and added to P 1 . Every time a new reaction is added, the maximum biomass production is computed given that all reactions in P 1 are knocked out. If the maximum biomass production is not equal to v ref biomass , the next reaction is drawn and added to P 1 . If the maximum biomass production is equal to v ref biomass , no more reactions are added to P 1 .

Step 3:

Step 2 is repeated n times to compute P 11 , P 12 , ..., P 1n . Then, the smallest (with regard to the number of reactions it contains) valid knockout set P 1x is chosen. If P 1x and P 0 are identical, no smaller knockout set can be chosen and the procedure is stopped. Otherwise, P 1x is set to be the new P 0 and Steps 2 and 3 are repeated.

It is not necessary to recompute the subsets 1 and 2. Subset 1 always remains the same and the new subset 2 will be a subset of the previous subset 2 containing only the reactions that are present in P 1x . Furthermore, it is important to observe that this approach is based on using the biomass production as main indicator because in our case study, the target production depends heavily on the obtainable biomass production. Therefore, it is not necessary to also compare the target production after each added reaction which saves computation time. In different cases, it might be advantageous to chose another main indicator or to compare both values after each newly drawn reaction The exploration is random and greedy and it is not guaranteed to obtain the smallest possible knockout set. The probability to achieve a small knockout set can be increased by choosing a high n When drawing reactions for the next knockout sets, if r 2 and r 3 are getting picked before r 1 , the limit for the biomass production might already be reached before r 1 is added because by knocking out r 2 and r 3 , the flux cannot deviate from the hyperpaths via r 4 . If this knockout set gets chosen for the next iteration, it is not possible anymore that r 1 gets picked anymore instead of r 2 and r 3 which would reduce the size of the knockout set.

and by repeating the whole exploration several times and comparing the results. Since the approach is greedy, always the smallest knockout set is chosen for the next iteration. It is likely that there are cases where choosing a different knockout set will lead to a better final result because it contains a better combination of reactions. An example is shown in Figure 2.4.

After applying this random exploration to the case study that will be presented in the results, it becomes clear that this described random exploration is too straightforward and it would be better to modify it in order to make it less greedy. The main idea is to not always take the smallest knockout set as next candidate but to introduce probabilities that guide the choice of the next candidate and aid in better exploring possible solutions. Still, smaller solutions should be prioritized.

Steps 1 and 2 of the modified version remain unchanged.

Step 3 is modified in the following way: Instead of repeating Step 2 n times, it is only done once to compute one new knockout set P 1 . All already computed knockout sets can be potential candidates for the next iteration. All computed knockout sets are unique. This means that if the newly computed knockout set was already computed previously, it is not added to the pool of all knockout sets. All computed knockout sets are sorted by size and split into buckets. In the following explanation, s is the size of the currently smallest knockout set. Furthermore, B is the size that is chosen for the buckets. This means that the first bucket b 0 contains all knockout sets with a size in between s and s + B -1 and the bucket b n contains all knockout sets with a size in between s + n • B and s + (n + 1) • B -1. However, there is a further adjustment for the first bucket. It contains all knockout sets of size s and additionally, the same amount of larger solutions. For example, if there are 50 knockout sets of size s, it also contains the 50 next bigger knockout sets. Furthermore, it should have a minimum size F . If there are less than F/2 knockout sets of size s, it will be filled with the next bigger knockout sets until it contains a total of F knockout sets. If in the current iteration, less than F knockout sets have been computed in total, the first bucket contains all knockout sets. Consequently, the limits for the subsequent buckets have to be adjusted. Given that the largest size that is still contained in bucket b 0 is s * , the bucket b n contains all knockout sets with a size in between (s * + 1) + n • B and

(s * + 1) + (n + 1) • B -1.
To choose the next candidate, first a bucket is chosen and afterwards a knockout set in the chosen bucket is picked. Starting with b 0 , each bucket has a probability p 1 to be accepted. If the bucket is rejected, the next bucket is tried. If a bucket is accepted, a knockout set from this bucket is drawn.

In this case, a uniform distribution is used which means that all knockout sets inside one bucket have the same probability.

After a knockout set is chosen, steps 2 and 3 are repeated. After each iteration, it is verified if a smaller knockout set was found. If after a specified amount of iterations, no smaller knockout set was computed, the random exploration stops.

The main idea behind the buckets is to make it easier that larger knockout sets get chosen as candidates but at the same time knockout sets of same or similar size have the same probability of getting picked. Additionally, since all computed knockout sets can be potentially drawn in the next iteration, this approach is less greedy than the first version where it was not possible to pick a previous candidate in a later iteration.

Implementation

Since in the second part, only single-objective optimization problems are solved, PolySCIP is no longer needed and all linear programs are solved with CPLEX 12.71 (IBM, 2016). They are implemented in C++ and the metabolic network was modeled using the library metnetlib which is available on https://gitlab.inria.fr/erable/kirikomics/metnetlib.

Results

To apply our approach, we used the production of ethanol in yeast since ethanol is an interesting target chemical due to its use as biofuel but it is also toxic to yeast.

All computations are done using the yeast 5.01 model [START_REF] Heavner | Yeast 5-an expanded reconstruction of the saccharomyces cerevisiae metabolic network[END_REF]. The network model contains a total of 2109 reactions and 2759 metabolites. After the removal of blocked reactions which are all reactions that can only have a zero flux for the given lower and upper bounds (see Section 1.2.4) and the subsequent removal of isolated metabolites, the model contains 1165 reactions and 914 metabolites.

We will first present the observations on the computed Pareto front. Afterwards, the results for the first approach are shortly summarized before showing the results for the second approach that is based on the idea of the hyperpaths.

Critical reactions

To identify reactions that might be involved in increasing the resistance of yeast against ethanol, the growth of single gene knockout yeast strains under ethanol stress was compared to the wild type growth under the same conditions. If the biomass production was decreased in the modified strain, it is assumed that the knocked out gene is implicated in developing a resistance against ethanol. Genes that are known to participate in growth mechanisms were not tested since knocking these genes out would likely lead to a reduced growth but not because they are influencing the resistance against ethanol.

Depending on the amount that the growth decreased compared to the wild type, the genes were separated into two groups. Group 1 contains genes whose knock out resulted in a biomass production between 50% and 70% of the wild type biomass production. Knocking out the genes in group 2 led to less than 50% of the wild type biomass production. Hence, genes in group 2 might be playing a more significant role in building a resistance against ethanol than genes in group 1.

Afterwards, the gene identifiers were used to link these genes to reactions in the network file. Not all genes can be associated with reactions. It is also possible to have several genes linked with the same reaction and one gene connected with several reactions in the network file.

The reactions that can be connected to the genes in groups 1 and 2 are the critical reactions that are used to compute the toxicity resistance score. A total of 61 critical reactions could be identified.

After removing blocked reactions, the list could be reduced to 41 remaining critical reactions. The reduced list can be found in the Supplementary Table A.3.

Computation of the Pareto front

The critical reactions were used to compute different tradeoffs between biomass production, ethanol production and the score measuring the potential resistance against ethanol. The multi-objective optimization problem was solved using PolySCIP version 2.0 [START_REF] Borndörfer | Polyscip[END_REF]. The critical reactions in group 1 have weight 1.0 when calculating the resistance score, reactions in group 2 have a weight of 2.0. Only the extreme points of the Pareto front were computed which allows to gain a broad overview of the whole Pareto front (see Figure 2.5). Some extreme points have very similar ethanol and biomass production values but differ considerably in their ethanol resistance score. In this case, a point with a higher ethanol resistance is preferable.

A list of computed extreme points and their identification number can be found in the Supplementary Only the extreme points were computed. The lower bound for the biomass production was set to 0.1 and the lower bound for ethanol production was set to 10. For a clearer representation, the third dimension, the toxicity resistance score, is displayed by a color scale. Red points have a low toxicity resistance score, blue points have a higher score. It is interesting to see that there are some points that have very similar ethanol and biomass production values but differ in the resistance score. It shows that it might be possible to gain a considerably better resistance against ethanol by choosing slightly worse production values for ethanol and biomass.

First approach

An important remark for the results for the first approach is that they were obtained using the original network model. The preprocessing step was introduced at a later point in time when the idea based on the identification of the hyperpaths was developed. As described in Section 2.3.2, the number of candidates for knockouts was reduced. After eliminating all reactions without a gene association and all transport and exchange reactions, 660 reactions remained in the candidate pool. Afterwards, for the remaining candidates, the coupled groups were identified and only one representative for each group was kept in the pool which reduced the candidate pool to 492 reactions. Subsequently, biomass and target production in the wild type were verified for single knockout of the remaining candidates. After discarding all candidates whose knockout resulted in a zero biomass or target production, 435 reactions were left in the candidate pool. These steps are the same for all tradeoffs and do not have to be repeated. The last step is shown on the example of the tradeoff 1 (as presented in Table 2.1). Having verified if a candidate reaction can have a zero flux in at least one tradeoff flux distribution, 65 more reactions could be removed from the pool. Finally, the remaining 370 candidates were sorted into the three groups. The first group was empty (must have non zero flux in the wild type and zero flux in the tradeoff). The second group contained 28 reactions (must have non zero flux in the wild type and can have zero flux in the tradeoff). Therefore, the third group included the last 342 reactions. Single knockouts were enumerated for all extreme points. These are however never sufficient to enforce the desired tradeoff flux. When we computed the double knockouts, we realized that enumerating all solutions is not practicable because there are only a few double knockouts that actually make the optimization problem infeasible. Consequently, almost all pairs of candidates were enumerated as double knockouts. After each iteration, for the newly found knockout set, a new constraint was added which meant that the time it took to enumerate the next solution slowly increased.

To reduce the number of enumerated pairs, we used the groups to prioritize certain candidates. We enumerated pairs that needed to include at least one of the 28 reactions from group 2 (i.e. they could appear together in a pair). This led to the computation of around 10000 double knockouts. Hence, each of the 28 reactions was paired with almost all of the other remaining 369 candidates (groups 2 and 3). This showed again that the main problem is that almost all pairs of candidates are feasible and do not disrupt the tradeoff flux. The underlying MILP is not restrictive enough. It might be faster just to try all pairs and verify the resulting biomass and target productions because it would avoid to solve the MILP for each new pair which renders this approach impracticable.

Second approach

In the following, the results for the second approach that is based on the hyperpaths are shown. Before applying the random exploration to identify knockout sets, the similarities and differences between the hyperpaths and tradeoffs are investigated.

Active reactions

In the subsequent computations of the hyperpaths, the number of active reactions was always limited to 500 reactions to prevent that too large hyperpaths can be enumerated.

To gain an overview of different flux distributions that are possible in the network, for 47 different tradeoffs, one hyperpath was computed and the active reactions were compared (see Figure 2.6). More than 300 active reactions were present in all hyperpaths. Since ethanol and biomass are produced in all of these 47 tradeoffs, reactions that are specific to these pathways have to participate in the solution. Although their values for biomass and ethanol production can differ significantly, topologically these solutions have many reactions in common. Furthermore, there are more than 100 reactions that appear in five or less hyperpaths. It is already an indicator that there might be several reactions that can be easily added or removed from a hyperpath without changing the resulting flux significantly. This impression is confirmed when comparing different hyperpaths for one tradeoff. The results for 1000 different hyperpaths for tradeoff 1 are shown in Figure 2.7. There are more than 200 reactions that are participating only in less than 10 hyperpaths. Hence, they do not seem to have a big impact on the solution and are probably only part of the hyperpath because it has to differ from previously computed hyperpaths that are then excluded by constraints. Besides, more than 350 reactions appear in all solutions.

Comparison of 10 hyperpaths for each of the 47 tradeoffs

Comparing the computed tradeoffs, some seem to be preferable because they have a higher ethanol production value and a better resistance score than others. However, it is necessary to compare many hyperpaths for each tradeoff because knocking out all outgoing reactions of a hyperpath does usually not enforce the exact values of the tradeoff. Since many reactions are essential for biomass production and ethanol production, the remaining subnetwork is still quite large and thereby also contains a lot of variability. In total, 884 reactions were active in at least one of these 1000 hyperpaths. Most of them participate either in all of the 1000 hyperpaths or they appear in just a very few cases. The binwidth for this histogram was set to five which means that, for example, the first bar counts all reactions that appeared exactly zero to four times in all 1000 hyperpaths. However, none of the represented reactions appeared zero times because reactions that were not part of any hyperpaths were not considered for this plot.

Before applying the random exploration approach to compute knockout sets, an adequate hyperpath has to be chosen. Thus, in a first step, for each tradeoff, 10 different hyperpaths were computed. All hyperpaths are evaluated by their ethanol production and resistance score when the biomass production is maximized as explained in 2.4.1.

First of all, it can be observed that ethanol production is dependent on the maximum biomass production. A high ethanol production can only be achieved if the maximum biomass production is low. If a high maximum biomass production is possible, it prevents a profitable ethanol production (see Figure 2.8). This happens because both biomass and product are directly dependent on the same carbon sources, which means that if all carbon molecules uptaken from the substrate are used for biomass production, there is nothing left to produce ethanol and vice versa. This fact can be exploited during the random exploration.

Moreover, different hyperpaths can lead to the same production values (see Table 2.2 for some selected results). For some tradeoffs, all ten hyperpaths have the same results (e.g. tradeoffs 9 and 11). In most of the other cases, the results for all ten hyperpaths are at least very similar (e.g. tradeoffs 2 and 4). Problematically, there are also instances where the results between the different hyperpaths can differ considerably. For example, the first nine hyperpaths for tradeoff 1 have very low maximum toxicity resistance scores and the achievable ethanol production is only mediocre. However, the tenth computed hyperpath differs in an important way and is superior in all values. Even the biomass production is slightly higher. This shows how important the choice of a specific hyperpath is. Enumerating ten hyperpaths for a tradeoff is not sufficient to ensure that favorable hyperpaths are available.

To show how the result is influenced if the toxicity resistance score is not taken into account, for each tradeoff, 10 hyperpaths were computed where only biomass production and target production were fixed in the optimization. The constraint (Equation 2.55) that is fixing the toxicity resistance to its optimal value from the tradeoff is removed from the optimization problem. Ideally, when including the toxicity resistance score in the optimization problem, the obtained hyperpaths should have a better toxicity resistance score than when it is omitted. Not fixing the toxicity score will not prevent the solver to choose a flux distribution that results in a high toxicity resistance score. However, it also does not encourage it and therefore, the probability for choosing a hyperpath with a lower toxicity should be higher.

The maximum toxicity resistance scores for both cases are displayed in Figures 2.9 and 2.10. Comparing the results, globally, higher maximum resistance can be achieved when the score is included in the optimization problem. Omitting it, only two of the resulting hyperpaths (tradeoff 19) have a maximum toxicity score over 20 as opposed to all ten when it is included. By including the score when computing the hyperpath, this can be achieved more consistently.

However, there are a few cases (e.g. tradeoff 4) where the ten hyperpaths computed without taking the toxicity resistance score into account have better maximum scores. A higher number of hyperpaths might be necessary to investigate if for this tradeoff a better toxicity resistance score can be reached when omitting the score from the optimization problem. It is possible that by enumerating more hyperpaths, other solutions are found where the score is at least as good as the highest one that can be obtained when the toxicity score is omitted from the optimization problem. Indeed, when enumerating 1000 hyperpaths for tradeoff 4, it is possible to obtain some hyperpaths that have a better toxicity score than any of the hyperpaths that are computed when omitting the score. For 47 tradeoffs with 10 different hyperpaths each, the resulting maximum biomass production is plotted against the maximum ethanol production that is possible. An ethanol production of more than 14 can only be achieved if the maximum biomass production is lower than 0.2.

Comparison of 1000 hyperpaths for chosen tradeoffs

Therefore, to get a more detailed view, 1000 different hyperpaths are enumerated for specific tradeoffs. Tradeoff 1 was chosen to investigate the fact that its tenth hyperpath differed significantly from the other nine and because it showed promising values for ethanol production. Additionally, the tenth hyperpath also had a high maximum toxicity score. Tradeoff 3 also had high values for ethanol production and toxicity score. One of the main objectives is to find a hyperpath with a very high ethanol production value and a high toxicity resistance score. Hence, both tradeoff 1 and tradeoff 3 needed to be explored further. Tradeoff 6 was chosen because the ten hyperpaths were inferior to the ones from tradeoffs 1 and 3 but we wished to explore whether enumerating a higher number might lead to better results. Tradeoff 13 was chosen for the same reason. Tradeoff 47 was selected because it is the tradeoff with the highest biomass production.

The results for tradeoff 1 are shown in Figure 2.11. First of all, by enumerating 1000 hyperpaths, it was possible to find more hyperpaths that have high ethanol production and high toxicity resistance score. Before, the tenth hyperpath was definitely the best result for tradeoff 1. It has a minimum target production of 17.10 and a maximum score of 21.01. The same result can be obtained for one of the newly enumerated hyperpaths. Additionally, several other hyperpaths were computed that have a slightly lower minimum ethanol production (in between 16.5 and 16.9) but also a slightly higher toxicity resistance score of around 22. Therefore, they are also interesting candidates for the random exploration.

Compared to tradeoff 1, the results for tradeoff 3 were inferior (see Figure 2.12). The highest minimum ethanol production is at 17.09. However, the corresponding maximum toxicity resistance score is only 13.01. The highest found maximum toxicity score was around 21.01 and the matching minimum ethanol production at 16.97. These values are slightly inferior to the best results of tradeoff 1.

For tradeoff 6, the 1000 enumerated hyperpaths did not contain any with a high maximum toxicity resistance score (see Figure 2.13). This might be also due to the fact that the optimal toxicity resistance score for tradeoff 6 was clearly lower than for tradeoff 1 or 3. Similarly, for tradeoff 13 which has a higher optimal toxicity resistance score than tradeoffs 1 and 3, it was possible to obtain hyperpaths with maximum toxicity resistance scores that were higher than any of the scores for tradeoffs 1 and 3. However, the corresponding minimum ethanol productions are lower (see Figure 2.14) because the optimal ethanol production of tradeoff 6 is also significantly lower than for tradeoff 1 or 3.

As expected, the computed hyperpaths for tradeoff 47 do not have very high values (see Figure 2.15). Especially, the ethanol productions are very limited. Interestingly, the minimum ethanol productions and the maximum toxicity resistance scores are not very diverse. It might show again that the variability of the flux distributions that remain possible is limited if the biomass production is high.

All in all, the results of these five different tradeoffs confirm that choosing a tradeoff with a high optimal toxicity resistance score can lead to hyperpaths with higher minimum toxicity resistance scores. Likewise, choosing a tradeoff with a lower optimal toxicity resistance score limits the minimum toxicity resistance scores for its hyperpaths. Furthermore, selecting a tradeoff with a high optimal ethanol production allows for the computation of hyperpaths with a high minimum ethanol production. Consequently, even though the computed hyperpaths have values for ethanol production and toxicity resistance score that differ from the optimal values of the corresponding tradeoff, the choice of the tradeoff does influence the result. It is important to choose a tradeoff whose optimal values are very close to the desired ones for target and biomass production and the toxicity score.

To justify the importance of the toxicity resistance score, the minimum ethanol production and maximum toxicity resistance score for hyperpaths including and omitting the toxicity resistance score in the optimization problem are compared. In tradeoffs 1 and 3 (see Figures 2.11 and 2.12), by including the toxicity resistance score, higher minimum toxicity resistance scores can be reached.

Discarding the toxicity resistance score, a slightly higher ethanol production can be achieved.

In the case of tradeoff 47, the toxicity resistance scores are better when ommitting the score (see Figure 2.15). However, as already mentioned before, fixing the toxicity score when enumerating the hyperpaths is restricting the model. Hence, all solutions that are feasible in this case are also feasible when omitting the toxicity score. It is therefore possible to obtain hyperpaths with a high toxicity score when the score is omitted.

The results for tradeoff 13 are similar to those for tradeoffs 1 and 3 but the difference between the highest minimum ethanol productions is larger. In contrast to tradeoffs 1 and 3, the results for tradeoff 6 are different because similar minimum toxicity resistance scores are obtained for both cases. This can be explained by the fact that tradeoff 3 has a lower optimal toxicity resistance score. Accordingly, it influences the resulting hyperpaths less and they are more similar to the hyperpaths generated without taking the toxicity resistance score into account.

Overall, analyzing the different hyperpaths, it is preferable to include the toxicity resistance score in the optimization problem to obtain hyperpaths with a higher maximum toxicity resistance score and similar minimum ethanol production.

Random exploration results

Afterwards, we used the random exploration presented in Section 2.4.2 to find a subset of the outgoing reactions that can be used as knockout set. To compare the values for the biomass productions that are computed during the random exploration to v ref biomass , a small tolerance of = 0.005 was used. This means that no new reactions are drawn for a knockout set if the maximum possible biomass production is smaller than or equal to v ref biomass + . The bin size was set to 10 and the minimum size for the first bin was set to 50. In each iteration, one new knockout set was computed. The enumeration was stopped if no new smaller knockout set was found after 1000 iterations. When choosing the candidate for the next iteration, the probability to accept a bin was set to 0.5 to increase the chance that not only the smallest current knockout sets are picked as candidates and that more diverse solutions are explored. Based on the results of the analyses of the different hyperpaths, knockout sets were enumerated for different hyperpaths of tradeoff 1 that were chosen based on their minimum ethanol productions and maximum toxicity scores (see Table 2.3). A total of 13 different hyperpaths were chosen and for each of these hyperpaths, the random exploration was repeated twice. After all runs were finished, interesting knockout sets were selected as follows. Only knockout sets that contained less than 25 reactions were selected. Furthermore, the minimum ethanol production needed to be above 14 and the maximum toxicity score above 20. For each run, only the smallest knockout sets were selected. However, if there was a slightly larger knockout set that led to an increase in the minimum target production, it was picked in addition. For example, if for one run the two smallest knockout sets have size 19 and both fulfill the desired threshold for ethanol production and toxicity score, they are selected. If there are knockout sets of size 20, they are not selected because the smaller knockout sets that are available are likely subsets of them. Therefore, the only reason to choose them is in the case where they lead to an increase in ethanol production and/or toxicity resistance score. In total, 18 knockout sets were selected in this way. The results are shown in Table 2.4. When selecting the smallest computed knockout sets, a first observation was that not all runs led to knockout sets that fulfilled the above mentioned criteria. In Table 2.4, we can see that in total, knockout sets got selected only from 12 out of the 26 runs that were done (two runs for 13 hyperpaths). Additionally, for several hyperpaths (ids 799, 800, 657, 31), neither of the two runs led to very small knockout sets. There are several possibilities to explain that. Since the proposed approach to select the knockout sets is randomized, it might be necessary to do more repetitions which might lead to smaller knockout sets also for these hyperpaths. Another possibility is that indeed, the outgoing reactions for these hyperpaths contain less favorable combinations of reactions and it is actually necessary to knock out more reactions to cut off the desired subnetwork. Lastly, the proposed parameter settings for the random exploration are the results of multiple test runs and seemed to offer a good tradeoff between time, exploring different (larger) solutions and driving the exploration towards smaller knockout sets. However, the possibility remains that further fine tuning of the parameters might improve the results. The 18 knockout sets contain 73 different reactions. A full version can be found in the Supplementary Table A.5. In Table 2.5, only the reactions are shown that occurred at least in five out of the 18 knockout sets. The three reactions that occurred the most (r 1110, r 0766 and r 1112) were always the reactions of subset 1 (reactions that have to be knocked out). The only exception was for hyperpath #169 for which both r 0766 and r 1112 were not part of subset 1. A possible explanation for this difference in relation to the other hyperpaths is that hyperpath #169 also had a slightly lower minimum ethanol production compared to all the other chosen hyperpaths (see Table 2.3).

The smallest knockout set that could be found and that satisfied the thresholds included 17 reactions. Problematically, many out of the 73 reactions represent transports in the organism (see Table 2.5) and it might be difficult to knock them out in vivo which means that the proposed knockout sets might not be practicable. As can be seen in Table 2.5, not all of the reactions do have associated genes which also makes it infeasible to knock them out in vivo. Furthermore, for example, reaction r 1110, which is an ADP/ATP transporter and was always part of subset 1 for the chosen hyperpaths and therefore played an important role in attaining the desired values for ethanol production and toxicity scores, has three associated genes. This means that it might require a greater effort to actually knockout this reaction in practice compared to other reactions that have less associated genes. The selected knockout sets were chosen based on the number of reactions that they contained. It might be necessary to select them based on the number of associated genes to identify knockout sets that need the smallest amount of gene knockouts. It is however not always necessary to knock out all associated genes of a reaction. If multiple associated genes for one reaction correspond to an enzymatic complex, knocking out one single gene might be enough to already disrupt its activity. For example, reaction r 1030 has five associated genes that are connected by 'and' in the model. Thus, knocking out one of the five genes should already interrupt the activity of this reaction.

Although certain reactions appeared in most of the selected knockout sets, a certain variability remained in the reactions that occurred less often which might make it possible to exchange at least some of the reactions that are more problematic or more costly (i.e. in terms of gene associations) to knock out. This also implies that by doing more runs or by broadening the selection of knockout sets, this choice could be further increased which could also help to identify knockout sets with fewer gene associations.

Discussion

We showed on the example of ethanol production in yeast that our developed approach can compute knockout sets that increase the target production and that ensure that reactions that are critical for a tolerance against the target can be active. To apply this approach, prior knowledge or experiments are essential to identify the reactions that are critical to improve the resistance of the microorganism against the toxic target. Without this information, the formulated multi-objective optimization problem is not applicable. If detailed knowledge for different reactions is on hand, their weights can be adjusted accordingly to improve the model. It is however also possible to modify the score that was used to capture the resistance against the toxic target if in specific cases a different model is preferable.

First of all, we computed the Pareto front to gain valuable insights about different efficient flux distributions in the network. As an advantage of this approach, we could see that there are some extreme points that have very similar ethanol and biomass production values but differ considerably in their toxicity resistance score. This already demonstrates that it is beneficial to include the toxicity resistance score in the optimization problem because desirable values for biomass and ethanol production can be reached without having reactions active that are critical for ethanol resistance.

Consequently, when focusing exclusively on biomass production and target production, it is possible that knockouts are proposed that will reduce the resistance against the toxic target. The Pareto front that displays tradeoffs between all three objectives can be used to choose more favorable outcomes for which knockouts should be enumerated.

Problematically, the MILP that was initially formulated to enumerate different knockout sets was not viable. It is possible to enumerate single knockouts but they are not sufficient to enforce the aspired flux values. For larger knockouts, the approach is unsuitable because the MILP is not restrictive enough. Moreover, even for single knockouts, the value of this approach is questionable because in retrospect, it seems to be easier to simply check all single knockouts instead of enumerating them.

As already proposed before, it might be possible to narrow down the enumerated knockouts to the ones that actually lead to a change by using a min-max formulation. The resulting increase in the complexity of the problem however might also make this approach less viable in practice, especially for larger networks.

In the second approach, we enumerated and analyzed hyperpaths for different points in the Pareto front. Interestingly, there was a huge overlap between all hyperpaths which confirmed that certain parts of the metabolic network are indispensable for the production of biomass and ethanol. Moreover, we could show the advantage of including the toxicity resistance score in our model. It must be said that by fixing the toxicity score when enumerating the hyperpaths, the model is more restricted. This means that all solutions that are feasible in this case are also feasible when omitting the constraint that fixes the toxicity score to a specific value. However, when omitting the score, the feasible space is larger and therefore, the solver is not guided to include the critical reactions in the hyperpath. Therefore, the probability that the computed hyperpaths lead to a good toxicity resistance score (based on the evaluation that was used) is lower. Analysing the different results revealed that multiple reactions that are critical for the resistance against ethanol are not essential for biomass and ethanol production.

One other benefit of the presented approach is that it is applicable to genome-scale networks. Although metabolic networks are commonly subject to analyses, many approaches are only applicable to smaller networks that represent just a condensed version of the metabolism and they struggle with larger networks due to the increased complexity. It is however important to analyse the complete metabolism to understand more complex relations and to model the organisms in a way that is as detailed as possible. Our method manages to extract knockout sets for very specific conditions from genome-scale metabolic networks.

The smallest computed knockout set contained 17 reactions. Realizing knockouts of this size in vivo is still challenging. However, we did not prove that 17 is the smallest knockout set that is possible.

The results that we obtained with the proposed parameter settings for the random exploration were promising but, of course, tuning the parameters further might actually also improve the results further, e.g. lead to the computation of smaller knockout sets. Additionally, we chose hyperpaths for the random exploration that aimed for a very high ethanol production. Giving more freedom and reducing these expectations might also allow for the computation of smaller knockout sets.

Up to now, the presented results are purely theoretical and their practical value remains to be evaluated biologically. Although the utilized metabolic network of yeast is very detailed and well curated, the simulated flux distributions for the introduced knockouts do not necessarily match reality. Furthermore, it is important to remark that it still needs to be investigated how applicable the proposed reactions are for knockouts. As mentioned during the presentation of the results, many of the reactions that were part of the smallest computed knockout sets were modeling different transports in the organism. It is known that these kinds of reactions are not always usable as candidates for knockouts in vivo. Some reactions also did not have associated genes. To avoid that reactions get included in the computed knockout sets that are unsuitable in practice, a list could be provided that contains only reactions that can be knocked out in vivo. This list can, for example, also be used to exclude reactions that have no known gene associations. By providing this type of list, the evaluation of the different hyperpaths also needs to be adapted because it will lead to outgoing reactions that cannot be knocked out. This means that the remaining variability of the hyperpath can be higher because it will be less restricted. Hence, different hyperpaths might actually get chosen as good candidates for the random exploration. One disadvantage of excluding certain reactions is that it might also bias the outcome and it is possible to miss important insights. So far, the approach models reactions knockouts. The smallest reaction knockouts do not necessarily correspond to the smallest gene knockouts. Hence, it might be preferable to adapt the approach to directly model gene knockouts.

It must also be noted that it is likely that certain genes have other functions that are not modeled in the metabolic network. Knocking them out might impact the results in unexpected ways. On the other hand, it could also mean that less knockouts are needed to achieve the theoretical results in vivo. It might therefore be advantageous to perform the knockouts in a sequential way. For example, knocking out one single gene from a reaction with multiple genes that correspond to an enzymatic complex might be enough to already disrupt the activity.

So far, we only considered knockouts in our approach. Other methods like OptForce [START_REF] Ranganathan | OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions[END_REF] and regulatory MCSs [START_REF] Mahadevan | Genome-scale strain designs based on regulatory minimal cut sets[END_REF] show that it is also interesting to identify up-and downregulation that will force a specific overproduction. In our case, the values that we wanted to enforce, which correspond to the different points in the Pareto front, were not obtainable in the enumerated hyperpaths. Indeed, when we maximized the biomass production and fixed it to its maximum for the consecutive minimization of ethanol production or maximization of the toxicity resistance score, the resulting values differed from the ones of the points in the Pareto front.

The subnetworks that remain after removing all of the outgoing reactions of a hyperpath include more than 300 reactions. This size also leads to high variability for the possible flux distributions. Especially the biomass production was often not restricted significantly which is problematic because we could also show that higher values for the ethanol production can only be obtained if the value for the biomass production does not exceed a certain limit. On the other hand, we also know that we cannot knockout any of the reactions on the hyperpath to further restrict the biomass production because the active reactions of the hyperpath were identified as vital to reach the values of the point in the Pareto front (i.e. for biomass production, ethanol production and toxicity resistance score). Since we cannot knockout any of the active reactions, one possibility could be to up-or downregulate some reactions on the hyperpath to reduce the undesired variability.

Another interesting aspect of this approach is that the overall framework is very flexible. We developed our approach to address situations where the desired target metabolite is toxic for the producing microorganism and might therefore result in a less efficient production. Hence, we computed tradeoffs between biomass production, ethanol production and a toxicity resistance score. It is however quite straightforward to modify this multi-objective optimization problem to account for different circumstances. For example, if the desired target metabolite is not toxic for the microorganism, the toxicity resistance score can be omitted from the problem. It would also be possible to compute tradeoffs between biomass production and the productions of two different target metabolites. Another example would be to consider a case where the production of the target metabolite also leads to the production of a toxic by-product. Here, the idea would be to identify tradeoffs that maximize the target production and minimize the production of the by-product. Independently of the adapted multi-objective optimization problem, hyperpaths can be enumerated accordingly, i.e. the newly formulated objective functions have to be fixed to their optimal values. Apart from that, the enumeration remains unchanged and also the concept for the random exploration does not change. Since our approach is suitable for genome-scale metabolic networks, this flexibility might make it interesting for different types of analyses.

Conclusion and perspectives

In this chapter, we presented a constraint-based approach that uses multi-objective optimization to identify tradeoffs between biomass production, target production and a score that measures the potential resistance of the microorganism against the toxic target. After tradeoffs have been computed, our method determines subnetworks that are needed to uphold the values of the three objectives.

The different subnetworks are evaluated based on the assumption that the biomass production is prioritized. The target production and the score are therefore always computed for the maximum possible biomass production. This assumption can be changed if a better model is on hand. Afterwards, our approach uses a random exploration to identify smaller knockout sets for promising subnetworks based on the outgoing reactions.

Applying our approach to the case-study of ethanol production in yeast, we were able to compute knockout sets with less than 20 reactions. It remains an open question to determine what the smallest possible size is. Additionally, we need more biological evaluation for the proposed knockout sets to establish how viable the suggested reactions are in practice.

The advantages of our approach are that it is applicable on genome-scale metabolic networks as we showed by using the yeast 5.01 model. Moreover, the framework is flexible and it should be possible to alter the objectives for different scenarios. Hence, we are also interested in applying our approach to other examples to further confirm its adaptability. We plan on polishing the approach and making it available as a tool on the Gitlab of the team. After obtaining some more biological opinion from collaborators on our results, we want to submit this work as a paper later this year.

Introduction

This chapter is based on the publication [START_REF] Ziska | Totoro: Identifying active reactions during the transient state for metabolic perturbations[END_REF] that was submitted to Bioinformatics.

The work is the result of a collaboration with Ricardo Andrade and Mariana Ferrarini. It is the continuation of an approach that Alice Julien-Laferrière developed during her PhD (Julien-Laferrière, 2016). This work was performed using the computing facilities of the CC LBBE/PRABI. The increasing availability of metabolomic data and their analysis are improving the understanding of cellular mechanisms and elucidating how biological systems respond to different perturbations [START_REF] Sevin | Biological insights through nontargeted metabolomics[END_REF]. Metabolomics can identify the metabolic capacities of an organism and this fact can be used to obtain a metabolic profile that characterizes the physiological response of a cell, tissue or organism to stress or to a general perturbation [START_REF] Roessner | What is metabolomics all about?[END_REF]. Different network-based strategies for metabolomic data analysis have been recently reviewed in (Perez de Souza et al., 2020) and amongst others, such strategies can be used to establish associations between metabolites or to integrate them into metabolic pathways. Metabolic profiles are often analyzed and interpreted with the help of bioinformatic software such as MetExplore [START_REF] Cottret | Metexplore: collaborative edition and exploration of metabolic networks[END_REF][START_REF] Frainay | Metaborank: network-based recommendation system to interpret and enrich metabolomics results[END_REF], MetaboAnalyst [START_REF] Xia | Metaboanalyst 3.0-making metabolomics more meaningful[END_REF][START_REF] Chong | MetaboAnalyst 4.0: towards more transparent and integrative metabolomics analysis[END_REF] or 3Omics [START_REF] Kuo | 3omics: a web-based systems biology tool for analysis, integration and visualization of human transcriptomic, proteomic and metabolomic data[END_REF]) that can identify the set of metabolites with a significant change in their concentration. The metabolomic data are projected on the annotated metabolic pathways in order to highlight the processes that may be linked to the observed changes.

The aforementioned software also try to integrate different kinds of omic data (such as transcriptomic, metabolomic or proteomic data) in order to give a deeper understanding of the studied mechanisms [START_REF] Cambiaghi | Analysis of metabolomic data: tools, current strategies and future challenges for omics data integration[END_REF]. Different approaches were reviewed in [START_REF] Rosato | From correlation to causation: analysis of metabolomics data using systems biology approaches[END_REF][START_REF] Ivanisevic | From samples to insights into metabolism: Uncovering biologically relevant information in lc-hrms metabolomics data[END_REF][START_REF] Stanstrup | The metarbolomics toolbox in bioconductor and beyond[END_REF] and software for the enrichment analysis of metabolomic data were evaluated and their results compared in [START_REF] Marco-Ramell | Evaluation and comparison of bioinformatic tools for the enrichment analysis of metabolomics data[END_REF]. However, metabolic pathways have subjective definitions and can differ between databases [START_REF] Ginsburg | Caveat emptor: limitations of the automated reconstruction of metabolic pathways in plasmodium[END_REF]. Additionally, this kind of analysis can make it hard to identify the connections between metabolites since they can be part of many pathways and it is thus possible to miss paths which traverse several pathways. Another approach is to use graph-based methods that allow us to consider the whole metabolism as an integrated system focusing on the parts that are connecting the metabolites of interest. Usually, these methods rely mainly on the network structure, chemical information and on an input list of metabolites [START_REF] Frainay | Computational methods to identify metabolic sub-networks based on metabolomic profiles[END_REF].

In (Acuña et al., 2012a;[START_REF] Milreu | Telling metabolic stories to explore metabolomics data: a case study on the yeast response to cadmium exposure[END_REF], a method is proposed in this direction that is based on the enumeration of metabolic stories. A metabolic story is defined there as the set of reactions that summarize the flow of matter from a set of source metabolites to a set of target metabolites and is characterized as a maximum directed acyclic sub-graph connecting the metabolites of interest. One of the drawbacks of this approach is that a metabolic story is acyclic and thus, it is not possible to obtain sets of reactions that contain cycles. However, cycles are common in metabolic networks and this assumption might thus not reflect reality. Additionally, it does not take the stoichiometry of the reactions into account. This can in turn lead to a set of reactions that is not feasible in practice.

Metabolite concentrations have been used to assess the responses to small perturbations in the context of constraint-based models [START_REF] Palsson | The challenges of in silico biology[END_REF][START_REF] Covert | Constraints-based models: regulation of gene expression reduces the steady-state solution space[END_REF][START_REF] Klamt | Stoichiometric and constraint-based analysis of biochemical reaction networks[END_REF]. In [START_REF] Reznik | Flux imbalance analysis and the sensitivity of cellular growth to changes in metabolite pools[END_REF], the authors used a method that is derived from the classical flux balance analysis (FBA) framework. They showed that the variables of the dual problem, the so-called shadow prices, which correspond to the sensitivity of FBA to imbalances in the flux, can indicate if a metabolite is a growth-limiting metabolite in FBA. In [START_REF] Rohwer | Identifying and characterising regulatory metabolites with generalised supply-demand analysis[END_REF][START_REF] Christensen | Tracing regulatory routes in metabolism using generalised supply-demand analysis[END_REF], methods are presented to identify regulatory metabolites and paths by varying in silico their known concentrations in a measured steady-state using supply-demand analysis. Therefore, these methods are based on the response of an organism to a relatively small perturbation and on the influence of the metabolite concentrations on the reaction rates of the system to return to the original equilibrium.

In this work, we focus not on the metabolite pools in one condition but on the difference of the obtained measurements between two conditions. We suppose that the difference of metabolite pools between two metabolic states can provide information on the transient state, that is, on the transition between the two measured conditions. In [START_REF] Sajitz-Hermstein | iremet-flux: constraint-based approach for integrating relative metabolite levels into a stoichiometric metabolic models[END_REF], the authors provide a method to integrate relative metabolomic measurements in order to make predictions about differential fluxes. They used a constraint-based approach which minimizes the distance between the two flux vectors of the two different states based on the ratio between the measured metabolite concentrations in both conditions. For both states, steady-state is assumed for the flux vectors. The authors identify differential fluxes between the two conditions whereas our approach will aim to find reactions that are potentially active during the transient state.

In [START_REF] Case | Reachability problems for continuous chemical reaction networks[END_REF], a similar problem was studied. The authors investigated reachability problems in chemical reaction networks. Given two different states of the network, the goal is to identify a path that leads the network from the first state to the second one. They prove that this problem can be solved in polynomial time. However, they also discuss that a variant of this problem in which the maximum size of the path is fixed is more difficult to solve. Our approach overcomes this limitation and is able to minimize the number of active reactions which is important since we are interested in identifying only the parts of the network that are potentially active during the transient state.

The method we propose uses a constraint-based modeling to enumerate sets of reactions that explain the changes in concentrations for some measured metabolites, i.e. how the system moved from a state to another. We implemented our approach in a software we called Totoro (for "Transient respOnse to meTabOlic pertuRbation inferred at the whole netwOrk level"), that is publicly available at https://gitlab.inria.fr/erable/totoro. It is implemented in C++ and depends on IBM CPLEX which is freely available for academic purposes. We also tested our method with data from pulse experiments with different carbon sources (glucose, pyruvate and succinate) in Escherichia coli.

Methods

A metabolic network can be represented as a directed hypergraph H(V, R, S) where V is the set of vertices, R the set of hyperarcs and S the stoichiometric matrix. Each c ∈ V represents a metabolite of the network and each hyperarc r ∈ R a reaction that connects two sets of disjoint metabolites Subs r , P rod r with Subs r , P rod r ⊆ V. The stoichiometric matrix S is a m × n matrix where each column represents a reaction and each row a different metabolite. It contains the stoichiometric coefficients which are positive if a metabolite is produced by a reaction and negative if it is consumed.

The set X ⊆ V contains all measured metabolites. The metabolomic data is given as a list which, for each measured metabolite in X, contains an interval. This interval describes by how much the internal metabolite concentration changed between two different states. Usually, small deviations for the measurements are available which can be used to calculate the minimum and the maximum possible difference between the internal metabolite concentrations. Furthermore, all reversible reactions of the network are split into forward and backward reactions.

We are interested in solving the following problem: Given a network H and a list containing the changes for some metabolite concentrations before and after a perturbation, we want to identify sets of reactions that were involved in diverting the system from the initial state before the perturbation to the state after the perturbation. We will present a constraint-based approach to solve this problem where the change of concentrations between two states is represented as an interval.

Core model

The variation of the concentrations in time of the metabolites in X can be written as:

dX dt = (S • v) X . (3.1)
In this equation, v is a flux vector and the (•) X operator means that only the entries of the vector corresponding to the metabolites in X are taken into account. We use [X] t to denote the concentration for the metabolites in X at time point t. Considering two points t 0 and t 1 in time and

∆ X = [X] t1 -[X] t0
, one can write:

∆ X = S • ϕ. (3.2)
In this case, each entry of the vector ϕ can be interpreted as the overall number of moles that passed through the reaction j during the time interval [t 0 , t f ] which corresponds to the area under the reaction rate curve in this time interval:

ϕ j = t1 t0 v j (t) • dt. (3.3)
Due to biological and technical variability that can arise from different replicates of the same experiment, we assume that the measured variations in concentrations of the metabolites in X are represented by an interval rather than using a fixed number:

∆ X = [∆ min X , ∆ max X ]. (3.4)
Furthermore, for the non-measured metabolites, we do not know if their concentration changed or not. Therefore, we will assume that a small variation is possible for all non-measured metabolites X = V \ X:

∆ X = [ min , max ]. (3.5)
Based on these assumptions, we can model the production or consumption of metabolites between two states by the following constraints:

∆ min ≤ S • ϕ ≤ ∆ max 0 ≤ ϕ j ≤ u j ∀j ∈ R. (3.6)
All ϕ j are positive and have an upper bound u j . ∆ min is a vector composed of ∆ min X and min while ∆ max is composed of ∆ max X and max .

The numerical values of the ϕ vector are difficult to interpret. The variable ϕ can only be zero or have a positive value. This means that we do not know if the activity of the corresponding reaction was increased or decreased during the shift compared to the baseline. We only know that if ϕ j is zero in the solution, reaction j was not active during the shift while if ϕ j has a non-zero value, reaction j was active during the shift. Hence, we are only interested in the reactions that have a non-zero ϕ because we want to identify the part of the metabolic network that was active during the metabolic shift. These reactions are represented by the support of the vector ϕ.

Minimizing the number of reactions and the variation of the concentrations for the non-measured metabolites

Since the number of possible paths that can explain the measured metabolic shifts can be very large, we will focus on finding the smallest solutions with regard to the number of active reactions that still explain the metabolic shift. This corresponds to the parsimonious assumption that the fewest possible resources are used or the smallest changes are made. Thus, we are interested in identifying minimum sets of reactions that play a major role in the metabolic shift. For each reaction j, a binary variable y j is then introduced that is set to zero if and only if the corresponding ϕ j is zero and therefore, the reaction is not part of the solution. In this way, these variables will correspond to the support vector of ϕ and it will be sufficient to minimize their sum:

y j = 0 ↔ ϕ j = 0 ∀j ∈ R y j ∈ {0, 1}. (3.7)
Additionally, to prevent that both a reaction j and its reversible j can be picked at the same time for one solution, the following constraint is used:

y j + yj ≤ 1 ∀(j, j) ∈ R. (3.8)
To minimize the number of reactions that are part of the solution, the objective function is written as:

min m j=1 y j . (3.9)
However, we are not only interested in minimizing the number of reactions in the solution but also in minimizing the variation in concentration for the non-measured metabolites X. Since the measured compounds are usually the more important ones for analyzing the biological experiment, it is reasonable to aim for solutions where other compounds do not accumulate or deplete a lot. This leads to the following minimization:

min |S • ϕ| X . (3.10)
On the other hand, we are trying to explain as much change in the concentration as possible for the measured metabolites:

max |S • ϕ| X . (3.11)
To combine both ideas in one objective function, a weight λ is used for both objectives:

min λ m j=1 y j + (1 -λ) |S • ϕ| X -(1 -λ) |S • ϕ| X . (3.12)
The value for λ should lie between 0 and 1. Finding a good balance between these two objectives can be challenging but necessary to identify meaningful biological solutions. This will be further discussed in Section 3.3.

Summing up, the mixed-integer linear program (MILP) that is implemented in our software Totoro is the following:

min ϕ,y λ m j=1 y j + (1 -λ) |S • ϕ| X -(1 -λ) |S • ϕ| X s.t ∆ min ≤ S • ϕ ≤ ∆ max 0 ≤ ϕ j ≤ u j ∀j ∈ R y j = 0 ↔ ϕ j = 0 ∀j ∈ R y j + yj ≤ 1 ∀(j, j) ∈ R
y j ∈ {0, 1}; λ ∈ (0, 1); u j , ϕ j ∈ R.

(3.13)

Enumerating different solutions

To enumerate different solutions, once a solution is found, it must be excluded for the next iteration. Two solutions are different if they do not contain the same reactions. We are using the following constraint where y * is a previously found solution vector: (3.14) This prevents that the exact same combination of reactions gets chosen again. Afterwards, we can solve the updated MILP again to compute a different solution. We repeat this process until no more new solutions can be found or until a desired number of solutions has been computed.

j∈R:y * j =1 y j ≤ m j=1 y * j -1.

Dealing with source/sink reactions and co-factors

In graph-based methods, it is known that looking for shortest paths without taking into consideration co-factors (for example ADP) can lead to irrelevant paths because such metabolites can introduce shortcuts through the network [START_REF] Frainay | Computational methods to identify metabolic sub-networks based on metabolomic profiles[END_REF]. In our case, although we are using a constraint-based approach and taking stoichiometry in account, similar problems can arise. When considering only shortest paths, depending on the presence of source or sink reactions and/or the value chosen for min and max , the active reactions in the solution can be highly disconnected. This makes them biologically less meaningful because it is not possible to identify possible pathways that played a role during the metabolic shift.

For example, if only the size of the solution is minimized, it is possible that changes in the concentration are just transferred to a close source or sink without actually selecting a pathway. To avoid this effect, it is important to block transport reactions in the network. Blocking transport reactions means that they cannot be part of a solution. However, if the substrates of a sink reaction are accumulated or the products of a source reaction are depleted in a solution, this indicates that the corresponding transport reaction can be part of the solution. Their use is limited by the chosen but it can be set to a very low or large value to imitate an infinite source or sink. Specific sources or sinks can be added to the problem by specifying a large negative ∆ min or a large positive ∆ max for certain metabolites. Therefore, transport reactions should always be blocked by setting their lower and upper bounds to zero.

A similar effect can happen if the value of epsilon for the non-measured metabolites is chosen too big. In this case, the changes in concentration of the measured metabolites can simply be distributed on (accumulated on or taken from) the nearby non-measured metabolites. This prevents that longer pathways are chosen which would actually connect several measured metabolites and could explain how the depletion of one measured metabolite leads to the accumulation of another measured metabolite (and vice-versa). However, this issue can be addressed by decreasing the value of λ in the objective function and thereby giving more weight to the function that minimizes the accumulation in non-measured metabolites. This should result in solutions that are larger but that connect the measured metabolites better than when only the number of reactions is minimized. Furthermore, it might be preferable to choose smaller epsilons to further restrict the accumulation/depletion of non-measured metabolites.

Before we minimized the accumulation and depletion of non-measured metabolites, we tried to prevent that shortcuts through the network are taken by limiting the amount of connected reactions that appear in the solution as active for each metabolite. Co-factors can usually have a very high degree which means that they are involved in many reactions. To prevent a high degree in the solution, the following constraint was used:

j∈R:i∈P rodj∪Subsj

y j ≤ D ∀i ∈ V. (3.15)
As a consequence, a metabolite cannot have more than D producing and consuming reactions that are part of the solution.

However, after including the minimization of the accumulation and depletion of non-measured metabolites in the objective function, we omitted this constraint. It was not necessary to treat co-factors differently since we can still obtain connected pathways as will be shown in Section 3.3.

Calculating the input deltas

To apply Totoro, it is necessary to provide the input deltas for measured metabolites. In most cases, they will not be available but they have to be calculated from the available metabolite concentrations of the different two conditions.

In the data that we used to apply Totoro, internal metabolite concentrations for the first condition (glucose baseline experiments) were given. They also included deviations that resulted from different replicates of the same experiment. The concentrations for the second condition (a pseudo-steady state after a pulse) had to be calculated from the metabolite concentrations of the first conditions and a fold change.

Therefore, for each measured metabolite x ∈ X, we computed the minimum (maximum) internal concentration for the pseudo-steady state [x] pss-min ([x] pss-max ) using the baseline glucose concentration [x] baseline with the given deviations d x and the fold changes f x :

[x] pss-min = ([x] baseline -d x ) • f x [x] pss-max = ([x] baseline + d x ) • f x .
(3.16)

Afterwards, we computed the minimum difference ∆ min

x and maximum difference ∆ max x between the internal concentrations of the glucose baseline and the pseudo-steady state:

∆ min x = [x] pss-min -([x] baseline + d x ) ∆ max x = [x] pss-max -([x] baseline -d x ).
(3.17) When the differences are negative, ∆ min x must have a distinct direction meaning that it is not possible that ∆ min x is negative and ∆ max x is positive. Therefore, measured metabolites with a fold change of 1.0 are considered as non-measured metabolites and are assigned the chosen generic . In some cases, for example if the fold change is close to 1.0 and the baseline deviations are large enough, it is possible that ∆ min x is negative and ∆ max x is positive. If one of them is clearly closer to zero, it will be set to zero. If no clear direction is determinable, the corresponding metabolite will be treated as a non-measured metabolite.

Results

To evaluate our approach, we used data from different pulse experiments with different carbon sources in E.coli as presented in [START_REF] Taymaz-Nikerel | Changes in substrate availability in escherichia coli lead to rapid metabolite, flux and growth rate responses[END_REF]. The authors measured the internal concentrations for several metabolites for a glucose baseline and for glucose, pyruvate and succinate pulse experiments. These data were used to apply the method on the E.coli core model [START_REF] Orth | Reconstruction and use of microbial metabolic networks: the core escherichia coli metabolic model as an educational guide[END_REF] and the E. coli iJO1366 model [START_REF] Orth | A comprehensive genome-scale reconstruction of escherichia coli metabolism-2011[END_REF] available from the BiGG database (King et al., 2015b). The E. coli core model consists of 72 metabolites and 95 reactions, the E. coli iJO1366 model of 1805 metabolites and 2583 reactions. We were interested in the difference between the glucose baseline and the pseudo-steady state shortly after the pulse experiment. In [START_REF] Taymaz-Nikerel | Changes in substrate availability in escherichia coli lead to rapid metabolite, flux and growth rate responses[END_REF], the authors provided the internal concentrations for the baseline, including the deviations for their measurements and the fold changes for the three different pseudo-steady states which we used to calculate the internal concentrations for each pseudo-steady state. In [START_REF] Taymaz-Nikerel | Changes in substrate availability in escherichia coli lead to rapid metabolite, flux and growth rate responses[END_REF], deviations for the measured concentration of the glucose baseline are given that were derived from several replicates of the same experiment. We used them to be able to calculate the minimum difference ∆ min X and maximum difference ∆ max X in the concentrations between the glucose baseline and each pseudo-steady state. A detailed explanation can be found in the Supplementary Material Section 3.2.5. The calculated ∆ min X and ∆ max X for all three pulse experiments can be found in the Supplementary Tables A. 6, A.7 and A.8. We used all measured metabolites that are present in the network and that had a significant change in their concentration as input. It should be noted that a change for each given metabolite must be either positive or negative. For further details, see the Supplementary Material Section 3.2.5. Furthermore, transport reactions cannot be chosen as part of the solution and therefore glucose, pyruvate and succinate were added as sources for the corresponding pulse experiment. Oxygen was added as another source because in [START_REF] Taymaz-Nikerel | Changes in substrate availability in escherichia coli lead to rapid metabolite, flux and growth rate responses[END_REF], the authors identified increased oxygen uptake rates during the pulse experiment. To allow unlimited growth, the biomass was added as sink.

The expected active reactions in the core metabolism of E. coli are displayed in Figure 3.1 for each pulse experiment.

E. coli core model

At first, the method was applied using the E. coli core model. To better understand how the different parts of our model impacted the solutions, we did several runs with different values for λ (0.1, 0.5 and 0.9) and (5 and 10) for each pulse experiment. Although a single solution should be enough to identify some pathways responsible for the shift, we wanted to see if we could also obtain alternative pathways. Furthermore, we wanted to investigate how the solutions evolve when they are slightly suboptimal. For each different parameter setting, 100 different solutions were therefore enumerated. The results are displayed using Escher (King et al., 2015a) in the Supplementary Figures A. 1 to A.15. In general, we could observe that solutions with λ = 0.1 are preferable since usually the goal is to have a final solution which is overall more connected. In this way, we were able to extract complete pathways that played a role during the metabolic shifts. This was the case for all three pulse experiments. A higher λ leads to solutions that are less connected since the optimizer prioritizes solutions with fewer active reactions in this case. However, this means that it is difficult to obtain complete pathways as solutions and it might be hard to interpret these solutions biologically. Nevertheless, the user is able to fine-tune the number of reactions in the final solution and the degree of connectivity (for instance, if the goal is to highlight only parts of the complete metabolic network instead of finding a connected pathway). We show this fine-tuning for the case of the glucose pulse, in which decreasing the parameter was used to obtain more connected pathways (see Figure 3.3).

By adjusting the parameters λ and , Totoro could propose complete pathways for all three pulse experiments. The predicted solutions did not use co-factors as shortcuts through the network. We therefore did not modify our method further to treat co-factors separately.

Pyruvate pulse

For the pyruvate pulse, we expected that the activity of the TCA cycle would go up and that reactions for gluconeogenesis would be active (see Figure 3.1). Both observations could be reproduced for λ = 0.1 for = 5 or 10. The TCA cycle was proposed to be active by Totoro in all the 100 enumerated solutions. The four measured metabolites citrate, isocitrate, L-malate and fumarate had positive input deltas and could thus be used as sinks. The results showed how the TCA cycle was fed from pyruvate either by the Phosphoenolpyruvate carboxylase (PPC) or by the combination of Pyruvate dehydrogenase (PDH) and Citrate synthase (CS). Furthermore, the pathway from pyruvate to D-fructose 6-phosphate was active, thereby also producing the biomass precursor glyceraldehyde 3-phosphate (G3P) in all of the 100 solutions. The pathway from pyruvate to G3P contains five reactions including the reversible reactions Enolase (ENO), Phosphoglycerate mutase (PGM), Phosphoglycerate kinase (PGK) and Glyceraldehyde-3-phosphate dehydrogenase (GAPD). Especially here, it is important to state that all these reversible reactions were predicted in the correct direction going from pyruvate towards G3P. One important difference between the results for = 5 and = 10 was that for = 5, the biomass reaction was chosen in all solutions which makes it slightly with their expected directions are highlighted in red whereas other non-essential reactions (but which nonetheless could be chosen) are depicted in grey. Each pulse is indicated by the short yellow arrow (Glc: glucose; Pyr: pyruvate and Suc: Succinate). During the glucose pulse, the glycolysis reactions (depicted in green) should be active in order to generate ATP from the hydrolysis of glucose. On the other hand, the pyruvate and succinate pulse experiments should show gluconeogenesis activation (also depicted in green but in the opposite sense), generating glucose-6-phosphate from these two carbon sources. Furthermore, the TCA cycle (depicted in blue) can be fed from pyruvate during the pyruvate and glucose pulses. During the succinate pulse, the overflow in the TCA cycle should lead to the production of pyruvate with a subsequent activation of gluconeogenesis to produce biomass precursors. The pentose phosphate pathway (depicted in purple) is most likely active in all pulses in order to generate biomass precursors; however, since this pathway is a mere interconversion of carbohydrates, there is no particular expectation as to the actual direction of these reactions. ). The metabolites that were given as input are highlighted in blue if the corresponding input deltas were below zero and red if they were above zero. Reactions that are highlighted in orange were chosen in almost all of the enumerated solutions. Reactions that are yellow were chosen only in around half of the solutions. White reactions were not chosen in any solution. The expected reactions of the gluconeogenesis and part of the TCA cycle are active in all 100 solutions. The reversible reactions of the gluconeogenesis are chosen in the correct direction. The figure was created using Escher (King et al., 2015a For λ = 0.9 (see Supplementary Figures A. 1 and A.2), neither the TCA nor the gluconeogenesis pathway were proposed to be active. Setting λ to 0.5 already improved the results: the TCA cycle was proposed as active but the gluconeogenesis pathway was only recovered in less than 50% of the solutions (see Supplementary Figures A.3 and A.4). We do not fix the objective value in our optimization problem after obtaining the first solution but in every iteration, the minimization problem is solved again after excluding the newly found solution. This means that the next solution can have the same objective value but it is also possible that the objective value is worse than in the previous iteration. In this particular case, the 100th solution had an objective value that is only 5.5% worse than the objective value of the first solution (see Table 3.1) which shows that, as concerns optimality, all 100 solutions were very similar. They also had very similar active reactions. Comparing the 100 enumerated solutions for λ = 0.1 and = 5, a total of 43 reactions with a specific direction were chosen in all solutions. Out of these 43 reactions, 24 were chosen in every solution (including reactions in the TCA cycle and the gluconeogenesis pathway). This means that certain core pathways were consistently picked also in slightly suboptimal solutions.

Looking at only the ten best solutions, already 38 out of the 43 reactions were identified. The missing reactions were mostly part of the pentose phosphate pathway which also contains reactions that were part of the solution only in a few cases. Even with only ten solutions, we were able to obtain the alternative pathways feeding the TCA cycle (PPC/PDH). This indicates that it is not necessary to enumerate a large amount of solutions to get significant results and to identify alternative pathways.

Glucose pulse

For the glucose pulse, we expected that reactions that are part of the glycolysis pathway would be active as they convert glucose into pyruvate generating energy. Consequently, the TCA cycle should also be fed (see Figure 3.1). For λ = 0.9 and 0.5, the active reactions proposed by Totoro were disconnected and it was not possible to identify active pathways. However, even for λ = 0.1 and = 5, only disconnected parts of the network were active (see Figure 3.3). Since we were interested in connected pathways, we decided to fine-tune the solutions by lowering the value of as much as possible. The result for = 1.2 can be found in Figure 3 the underlying optimization problem infeasible. As previously mentioned, decreasing the value of allowed us to obtain more connected solutions. For = 1.2, we got solutions that linked intermediate metabolites of the glycolysis pathway to the TCA cycle through the PPC reaction. Followed by a reversed Malate dehydrogenase (MDH), reversed Fumarase (FUM) and Fumarate reductase (FRD7), the solutions allowed for an accumulation of the input metabolites L-malate, fumarate and succinate. However, in some solutions, the TCA cycle was additionally fed by PDH and Citrate synthase (CS) to account for the accumulation of citrate. This means that when the solutions are disconnected and this is unwanted, better results might be obtained by lowering the value of parameter . Again, the 100 solutions were very similar (λ = 0.1, = 1.2). They accounted for a total of 47 reactions (with distinct directions) and 30 of these appeared in all solutions. Similarly to the pyruvate pulse, the difference in these solutions were mostly based on a few reactions that are not part of the main pathways (glycolysis/TCA cycle). One critical observation is that the D-glucose transport reaction (GLCpts) was not part of every solution although glucose should be used as important source. However, the bacteria were already grown in glucose as the baseline, which in turn might be a reason why glucose was already internalized prior to the initial pulse. When comparing the objective values for these 100 solutions, the absolute difference between the first solution and the 100th one was similar to the one observed for the pyruvate pulse (see Table 3.1). However, proportionally this value was 21.8% worse than for the first solution. When we repeated the run for λ = 0.1 and = 1.2 with 50 iterations, the D-glucose transport reaction was part of 42 solutions. For ten iterations, this reaction was picked in all ten solutions. Hence, the glucose transport reaction was active in solutions with the best objective values. This showed that although the solutions remained very similar, there was a decline in their quality. For the pyruvate pulse, we saw that it is not necessary to enumerate a large amount of solutions.

Succinate pulse

After the succinate pulse, part of the TCA cycle should always be active. Furthermore, the gluconeogenesis pathway should be active to produce G3P and glucose-6-phosphate from succinate. Again, the results for λ = 0.5 and 0.9 led to smaller solutions that were more disconnected (see Supplementary . Therefore, we focused on the analysis of the results for λ = 0.1 (see Supplementary Figure A.16 and Figure 3.4). For both = 5 and 10, succinate entered the TCA cycle and turned into oxaloacetate. Totoro proposed two possibilities to output the excess of the TCA cycle: Either phosphoenolpyruvate (PEP) was produced by PEP carboxykinase (PPCK) or by PEP synthase (PPS) using pyruvate as intermediate substrate. Subsequently, PEP was, as expected, transformed to G3P. The lower right part of the TCA cycle predicted as active can be explained by the fact that the concentration of L-glutamate decreased and the concentration of citrate increased. The active reaction in this part connected these two metabolites. Furthermore, reactions of the pentose phosphate pathway were proposed as active and the biomass precursors R5P, E4P and G3P were produced.

The results for = 5 and 10 were very similar. For example, one difference was that for = 10, the reverse D-lactate dehydrogenase (LDH) was predicted to be active in 56 solutions which led to a small accumulation of D-lactate. It does make sense biologically because in general, D-lactate is one of the main products of the fermentation but we do not have the measurements for the concentration of D-lactate for this pulse experiment to actually verify this observation. However, in total, the differences were negligible and in contrast to the glucose pulse, the parameter had a lower impact on the outcome. Again, the core reactions of all 100 solutions were very similar. In total, 41 reactions (with distinct directions) appeared in all 100 solutions (for λ = 0.1, = 5). We observed that 22 of these were always active (mostly in the gluconeogenesis pathway and part of the TCA cycle). The objective values for all 100 solutions were extremely close (see Table 3.1).

E. coli iJO1366 model

Based on the results for the E. coli core model, we only did runs with λ = 0.1 for the E. coli iJO1366 model. The inputs were updated because this network contains more metabolites and therefore, more measured metabolites could be added. The amount of iterations was decreased to ten because the runtime in the larger network is significantly higher and we had already established in the core model that it is not necessary to enumerate a larger amount of solutions. To decrease the runtime for each solution, CPLEX was configured differently. The relative MIP gap tolerance was set to 0.05 which means that the solver will stop an iteration if a solution is found that is within 5% of the optimal. This allows for a faster result and we could see in the core model that the first 100 solutions tend to be very similar. This means that even if we are enumerating slightly suboptimal solutions, we should be able to compute solutions that are very similar to the actual optimal solution. If the 5% limit is not reached after 48h, the iteration is stopped. The memory usage of CPLEX was limited to 10 GB. The runtime for the different pulse experiments differed a lot. The results for the pyruvate and glucose pulse were computed on a cluster. For the pyruvate pulse, the 5% limit was reached only in three iterations (see Supplementary Table 3.2). All other iterations were stopped after 48h. However, all solutions obtained were within 7% of the optimum. Thus, we still took them into account when analyzing the predicted active reactions. In none of the iterations for the glucose pulse, the 5% limit was reached. The obtained solutions were within 8.5% of the optimal value (see Supplementary Table 3.3).

In contrast to the pyruvate and the glucose pulse, the 5% limit was reached in all iterations for the succinate pulse and computing all ten solutions took less than 5 minutes on a personal machine (2.90GHz Intel i7-7820HQ CPU, 16GB RAM). This shows that the constraints describing the input deltas in the MILP have a large influence on the difficulty of the optimization problem, and thus also on the runtime. However, although the obtained solutions were suboptimal, the active reactions predicted by Totoro for the core metabolism are very similar to the best results of the E. coli core model for all three pulse experiments. The additional measurements that were added as input deltas for the large network were mostly amino acids (see Supplementary Tables A.6 to A.8). In [START_REF] Waschina | Metabolic network architecture and carbon source determine metabolite production costs[END_REF], the authors show for the example of amino acid production in E. coli how the production cost for individual amino acids can depend on the available carbon source, and reactions close to the entry point of the carbon source might have considerably higher fluxes (see Figure 3.5). Indeed, from the experimental data, valine only accumulated during the pyruvate pulse, and was depleted with the other two carbon sources. Pyruvate is a direct precursor for valine production. Therefore, we expected that reactions of the valine metabolism should play a greater role in the predicted results for pyruvate compared to the other two pulses. Totoro predicted a higher turnover from pyruvate to valine, which resulted in the accumulation of this amino acid. Even though the pathway was also predicted as active for the glucose pulse, it was consumed more in this case (see Supplementary Table 6). In accordance The solver stopped either if a solution within 5% of the optimal value was found or after 48 hours.

Only in three iterations, the 5% limit was reached. In all other iterations, the solver stopped after 48 hours. However, the obtained solutions had objective values with 7% of the optimum. Thus, we still took them into account when analyzing the predicted active reactions.

Iteration Objective value Optimality (%) Runtime (h) Results for glucose (λ = 0.1, = 1.0). In total, 10 solutions were computed. The table shows the objective value for each solution and how close this value is to the optimum (in %). The solver stopped either if a solution within 5% of the optimal value was found or after 48 hours. The 5% limit was never reached and the solver always stopped after 48 hours. All computed solutions were within 8.50% of the optimal value.

with the predictions in [START_REF] Waschina | Metabolic network architecture and carbon source determine metabolite production costs[END_REF], another example is the accumulation of threonine during the succinate pulse. As shown in Figure 3.5, threonine and succinate are closely connected, and

Totoro predicted active reactions leading to its accumulation during the succinate pulse. Compared to the results for succinate, Totoro predicted more active reactions consuming threonine during the glucose pulse and less active reactions producing threonine during the pyruvate pulse, resulting in the depletion of this amino acid in these two cases (see Supplementary Table 6). Moreover, only during the glucose pulse, phenylalanine was accumulated. As a result, Totoro proposed more reactions of the phenylalanine metabolism as active compared to the pyruvate and succinate pulses, in accordance with the predictions in [START_REF] Waschina | Metabolic network architecture and carbon source determine metabolite production costs[END_REF] of lower cost to produce this amino acid with glucose as carbon source (see Figure 3.5).

Another interesting pathway we noticed as active only in the glucose results was the murein recycling pathway. Murein (or peptidoglycan) is a polymer consisting of sugars and aminoacids and is a major component of cell wall in bacteria. As bacteria grow and the cells divide, this layer of peptidoglycans breaks and the fragments are transported back for reutilization [START_REF] Goodell | Recycling of murein by escherichia coli[END_REF][START_REF] Van Heijenoort | Peptidoglycan hydrolases of escherichia coli[END_REF]. Indeed, around 60% of the murein sacculus is thought to be cleaved and reused at each generation in E. coli [START_REF] Park | How bacteria consume their own exoskeletons (turnover and recycling of cell wall peptidoglycan)[END_REF]. Although merely speculative at this point, as glucose is the carbon source with the highest growth rate for E. coli, it is plausible to assume that the higher growth rate when compared to pyruvate or succinate results in a higher amount of murein to be recycled at each generation.

Discussion

Totoro was able to predict expected pathways as active based on the differences in the measured concentrations for some internal metabolites. We could show that in general, it is preferable to use λ = 0.1 though the method is not critically sensible to this setup, being robust to small perturbations. However, it is worth noting that a higher λ can lead to smaller solutions but they are biologically irrelevant. We were interested in extracting complete pathways that explain the changes in concentration between two different conditions. We could see that a reduction of could be used to obtain a more connected solution and therefore, in our case, a smaller led to better solutions. However, there might be situations where we are more interested in only local changes around the measurements. In this context, it might be advantageous to choose a higher . Furthermore, since we are doing more than just minimizing the number of active reactions, something which would be closer to looking for shortest paths, we did not encounter problems specific to co-factors. By splitting reversible reactions, Totoro is able to predict distinct directions for them. Both in the core network and in the larger network, we were able to recover pathways that make sense biologically. Additionally, although the larger network contains more reactions and we added more input deltas, the predictions for the core metabolism of E. coli were similar to the results for the core network. It must be however noted that the predictions do depend on the measured metabolites. If for large parts of the network, no metabolite concentrations are measured, Totoro will likely not be able to find active pathways for these parts of the network. Moreover, we could also see that it is not necessary to enumerate a high number of solutions which is especially important when larger networks are used and the runtime of Tororo increases. We enumerated 100 different solutions for the core network. However, in our case, the enumerated solutions are very similar and a large amount of reactions appeared in all 100 solutions. Therefore, already one solution would have been sufficient to infer the most important reactions that were proposed to be active.

Conclusion

In this work, we presented Totoro, a method that identifies active reactions during the transient state based on the differences in the concentrations for some measured metabolites in two different states and we showed its prediction power on the example of different pulse experiments in E. coli.

Our method Totoro only uses metabolomic data as basis for the prediction. It is able to handle full networks which take into account in the model stoichiometry, cycles, reversible reactions as well as co-factors. With the current developments, it gets more common to have different kinds of data available which creates a need for methods that combine, for instance, metabolomic, transcriptomic and proteomic data. In [START_REF] Pandey | Enhanced flux prediction by integrating relative expression and relative metabolite abundance into thermodynamically consistent metabolic models[END_REF], the authors proposed a framework that predicts differential fluxes. It consists of three different methods that use either thermodynamic and transcriptomic data or thermodynamic and metabolomic data or thermodynamic, transcriptomic and metabolomic data at the same time. It might then be interesting in the future to adapt our approach to be able to integrate different kinds of data.

Conclusions and perspectives

In this thesis, I presented two different methods for the analyses of metabolisms of microorganisms. One approach focused on the identification of knockout sets, the other one on the analysis of metabolic shifts.

The first method that I described was the main work of my PhD. It can be used to identify knockout sets that increase the production of a valuable target metabolite in a microorganism in the scenario where the target metabolite is toxic for the microorganism and its accumulation can therefore inhibit growth or lead to a decline in its production. In the given approach, the resistance of the microorganism against the toxic target was measured based on the activity of some critical reactions that were identified experimentally beforehand.

In the studied example of ethanol production in yeast, already in the first part of the approach where tradeoffs between biomass and ethanol productions and the toxicity resistance score are calculated, we could indeed see that there are flux distributions that have very similar biomass and ethanol productions but that differed significantly in the calculated toxicity resistance score. This showed that not all of the critical reactions are necessary for growth and ethanol production. Hence, accounting specifically for the toxicity is important to ensure that the microorganism can keep its natural resistance. We could further demonstrate the advantage of the toxicity resistance score by comparing different hyperpaths in the second part of our approach.

Applying the random exploration, we were able to identify smaller subsets of the outgoing reactions for some hyperpaths that can be used as knockout sets. However, we still need to obtain a more biological examination of the computed knockout sets to establish how viable they would be in practice.

One of the main advantages of our approach is that it is applicable to genome-scale metabolic networks. Moreover, the described framework is flexible and it should be possible to adapt it to identify knockout sets for different examples, not just cases where the target is toxic for the microorganism. We plan to submit our method as paper after receiving a biological view of our results from our collaborator. Furthermore, the implementation in C++ that uses PolyScip as solver for the multi-objective optimization problems and CPLEX for all other MILPs will be made available on the Gitlab of our group.

The second method that I presented in this thesis is called Totoro and it can be used to analyze metabolic shifts. Totoro predicts reactions that are active during the transient state that occurs after a perturbation. It minimizes the change in concentrations for unmeasured metabolites and also the number of active reactions during the transient state to account for a parsimonious assumption. It predicts distinct directions for reversible reactions. An implementation of Totoro in C++ that uses CPLEX as solver for the underlying MILPs is available at https://gitlab.inria.fr/erable/ totoro.

On the example of three different pulse experiments in E. coli, we could show that this constraintbased approach is also applicable to larger network models. We could reproduce the main observations obtained in the E. coli core model in the E. coli iJO1366 model. Furthermore, we showed that by limiting the accumulation/depletion of unmeasured metabolites in the MILP and also prioritizing their minimization in the objective function, Totoro can predict connected pathways. We did not encounter problems specific to co-factors.

With the current developments, it gets more common to have different kinds of data available. et al., 2015a).

TITRE en français

Modèles et algorithmes pour étudier et exploiter le métabolisme des micro-organismes RESUME en français Dans cette thèse, j'ai travaillé sur deux approches différentes pour étudier le métabolisme des microorganismes. La première méthode identifie les knock-outs qui augmentent la production d'une métabolite cible dans un micro-organisme dans le cas où le métabolite cible est toxique pour le micro-organisme utilisé. Dans la première partie de l'approche, un problème d'optimisation multiobjectifs est formulé qui calcule des compromis entre la production de biomasse, la production du métabolite cible et un score qui mesure la résistance possible à la toxicité du métabolite cible. Dans la deuxième partie, des knock-outs sont calculés sur la base de l'identification et de la séparation des sous-réseaux qui peuvent atteindre les valeurs de production souhaitées identifiées dans la première partie. L'approche est applicable aux réseaux métaboliques à l'échelle du génome, comme est montré dans l'étude de cas sur la production d'éthanol dans la levure. La deuxième méthode, appelée Totoro, a été développée pour l'analyse des changements métaboliques. Elle intègre les concentrations internes de métabolites qui ont été mesurées avant et après une perturbation dans les réseaux métaboliques. Il prédit les réactions qui étaient actives pendant l'état transitoire qui s'est produit après la perturbation. Totoro est une approche basée sur les contraintes qui prend en compte la stoechiométrie du réseau. La méthode est appliquée à trois expériences d'impulsions dans Escherichia coli pour montrer qu'elle peut récupérer des voies actives connectées et prédire des directions distinctes pour des réactions réversibles qui sont conformes aux observations biologiques connues. Totoro est applicable aux réseaux métaboliques à l'échelle du génome.

MOTS-CLEFS en français métabolisme; modélisation des réseaux métaboliques; énumeration; knock-outs; changement métabolique; hypergraphes dirigés; programmation sous constraintes Title in english Models and algorithms for investigating and exploiting the metabolism of microorganisms Abstract in english In this thesis, I worked on two different constraint-based approaches to study the metabolism of microorganisms. The first method identifies knockouts that increase the production of a target chemical in a microorganism in the scenario where the produced target metabolite is toxic for the utilized microorganism. In the first part of the approach, a multi-objective optimization problem is formulated that computes tradeoffs between biomass production, target production and a score that measures the possible toxicity resistance against the toxic target metabolite. In the second part, promising knockout sets are computed based on identifying and cutting off subnetworks that can lead to the desired production values identified in the first part. The approach is applicable to genome-scale metabolic networks which is shown in the case-study of ethanol production in yeast.

The second method which is called Totoro was developed for the analysis of metabolic shifts. It integrates internal metabolite concentrations that were measured before and after a perturbation into genome-scale metabolic networks. It predicts reactions that were active during the transient state that occurred after the perturbation. Totoro is a constraint-based approach that takes the stoichiometry of the network into account. The method is applied to three pulse experiments in Escherichia coli to show that it can retrieve connected active pathways and predict distinct directions for reversible reactions that are in accordance with known biological observations. Totoro is applicable to genome-scale metabolic networks.

Keywords in english metabolism; metabolic network modeling; enumeration; knockouts; metabolic shift; directed hypergraphs; constraint-based programming

Contents 3 . 2 . 5

 325 Calculating the input deltas . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.3.1 E. coli core model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.3.2 E. coli iJO1366 model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Objective space of a multi-objective linear program for k = 2. In this small example, the set of points {a, b, c, d, e, f, g, h} should represent the full set of feasible points Z. (a) All nondominated points Z N = {a, b, c, d, e, f } are highlighted in green. The dotted lines illustrate that it is not possible to improve one of the objectives without decreasing the other objective. The points {g, h} are dominated. For example, compared to point c, the objective values in point h are both smaller and therefore, h is dominated by c. The set of all nondominated points is also called Pareto front. (b) The nondominated points can be divided into supported (highlighted in green) and unsupported (highlighted in red) nondominated points. Point b is unsupported because it is dominated by an infeasible convex combination of a and c which is illustrated by the dotted line.

Figure 1 . 3 :

 13 Figure 1.3: Different compound graphs. The figures were adapted from[START_REF] Lacroix | An introduction to metabolic networks and their structural analysis[END_REF]. (a) An undirected graph. The edges represent reactions that transform one compound into another one. It is difficult to represent reactions that have multiple substrates and/or products. This means that the three displayed edges A ↔ B, B ↔ C and B ↔ C could correspond to three different reactions but they could also be the result of, for example, the reaction B ↔ A + C + D. Since no distinct directions can be represented, it is not possible to distinguish between reversible and irreversible reactions. (b) An undirected hypergraph. A hyperedge has one arbitrary subset of vertices. Therefore, it is not possible to distinguish between the substrates and products of a reaction. Here, the represented reaction could be B ↔ C + D but other interpretations like C ↔ B + D are also possible. Additionally, it is not possible to represent that a reaction is irreversible. (c) A directed hypergraph. The possible directions for each reaction are displayed by arrows. This representation allows for a simple discrimination between reversible and irreversible reactions. Additionally, there is a clear separation between substrates and products. Reaction A ↔ B is reversible whereas B → C +D is irreversible.

Figure 1 . 4 :

 14 Figure 1.4: Directed hypergraph and stoichiometric matrix representation. (a) A metabolic network can be represented as a directed hypergraph. In a directed hypergraph, the set of reactions R = {r 1 , r 2 , r 3 , r 4 , r 5 , r 6 , r 7 } is represented by the hyperarcs. The set of metabolites M = {S 1 , S 2 , A, C, B, T } is represented by the vertices. It is possible to display the stoichiometric values for each reactions by labeling the hyperarcs. In the given example, all reactions except reaction R 3 are irreversible. (b)The network can also be represented as stoichiometric matrix S. The rows of S correspond to the metabolites in M and the columns to the reactions in R. The stoichiometric values are displayed as coefficients of the matrix. If a coefficient is positive, the corresponding metabolite is produced by the reactions. If it is negative, it is consumed. Hence, in the matrix representation, a direction is assumed for each reaction. However, it is not possible to see if a reaction is reversible or not. Lower and upper bounds for reactions have to be stored separately.

Figure 1 . 5 :

 15 Figure 1.5: Flux balance analysis. The figure was taken from[START_REF] Kim | Methods for integration of transcriptomic data in genome-scale metabolic models[END_REF]. (a) Example of a small metabolic network. Reaction 1 is a source reaction that imports metabolites into the boundaries of the system. Reaction 3 exports metabolites. The network can be represented as stoichiometric matrix. Since reaction 1 can uptake metabolites from an external source, is does not have a substrate. In the same way, reaction 3 does not have any products. (b) Often, steady-state is assumed which means that concentration for metabolites does not change. (c) The feasible flux space can be restricted by adding lower and upper bounds for reactions.

Figure 1 . 6 :

 16 Figure 1.6: Different types of flux coupling. The figure was taken from (Burgard et al., 2004). Different scenarios for flux coupling are shown based on the flux ratio limits of R min = min v 1 /v 2 and R max = max v 1 /v 2 which have to be computed.

9 Figure 2 . 1 :

 921 Figure 2.1:Example for a hyperpath. All reactions that have a non-zero flux in the solution and participating metabolites are highlighted in green. The active reactions are r 9 , r 4 , r 3 , r 10 . The only internal metabolite is T arget. The metabolites S 1 , S 2 , C, D, E, Biomass are on the border of the hyperpath. Reactions r 2 , r 8 , r 9 are incoming and reactions r 1 , r 5 , r 9 are outgoing. Given that reaction r 3 is a critical reaction, using reactions r 9 , r 4 , r 3 instead of r 1 , r 2 is increasing the resistance against the toxic target.

6 Figure 2 . 3 :

 623 Figure 2.3: Topological precursors. The figure was adapted from(Acuña et al., 2012b). The metabolites S 1 , S 2 , S 3 and S 4 are the sources of the metabolic network. Metabolite T is the target node whose topological sources should be identified. For example, if M 0 = {S 1 , S 2 , S 3 }, it is possible to produce A because both substrates of reaction r 1 are present. Therefore, M 1 = {S 1 , S 2 , S 3 , A}. Now, both substrates of reaction r 2 are available and C can be produced which means that subsequently also T can be produced. Afterwards, no more changes can be made and a fixed point is reached. Hence, F wd({S 1 , S 2 , S 3 }) = {S 1 , S 2 , S 3 , A, C, T }.

4 Figure 2 . 4 :

 424 Figure 2.4: Disadvantage of the first version of the random exploration. In this example, reactions r 1 , r 2 and r 3 are outgoing reactions and the metabolites A, B and C are metabolites that lie on the border of the hyperpath. The metabolites D, E and F and reaction r 4 are outside of the hyperpath.When drawing reactions for the next knockout sets, if r 2 and r 3 are getting picked before r 1 , the limit for the biomass production might already be reached before r 1 is added because by knocking out r 2 and r 3 , the flux cannot deviate from the hyperpaths via r 4 . If this knockout set gets chosen for the next iteration, it is not possible anymore that r 1 gets picked anymore instead of r 2 and r 3 which would reduce the size of the knockout set.

Figure 2 . 5 :

 25 Figure 2.5: Computed tradeoffs between biomass production, target production and toxicity resistance score. Only the extreme points were computed. The lower bound for the biomass production was set to 0.1 and the lower bound for ethanol production was set to 10. For a clearer representation, the third dimension, the toxicity resistance score, is displayed by a color scale. Red points have a low toxicity resistance score, blue points have a higher score. It is interesting to see that there are some points that have very similar ethanol and biomass production values but differ in the resistance score. It shows that it might be possible to gain a considerably better resistance against ethanol by choosing slightly worse production values for ethanol and biomass.

Figure 2 . 6 :

 26 Figure 2.6: Histogram of active reactions in 47 different tradeoffs. For 47 extreme points, one hyperpath was computed for each. The figure shows in how many hyperpaths a reaction was active among these 47. A large part of the active reactions was present in all of them.

Figure 2 . 7 :

 27 Figure 2.7: Histogram of active reactions in 1000 different hyperpaths for tradeoff 1.In total, 884 reactions were active in at least one of these 1000 hyperpaths. Most of them participate either in all of the 1000 hyperpaths or they appear in just a very few cases. The binwidth for this histogram was set to five which means that, for example, the first bar counts all reactions that appeared exactly zero to four times in all 1000 hyperpaths. However, none of the represented reactions appeared zero times because reactions that were not part of any hyperpaths were not considered for this plot.

Figure 2 . 8 :

 28 Figure 2.8: Dependency of the maximum ethanol production on the maximum biomass production.For 47 tradeoffs with 10 different hyperpaths each, the resulting maximum biomass production is plotted against the maximum ethanol production that is possible. An ethanol production of more than 14 can only be achieved if the maximum biomass production is lower than 0.2.

Figure 2 . 9 :

 29 Figure 2.9: Maximum toxicity resistance score for hyperpaths where the score was included in the optimization problem.For each of the 47 tradeoffs, the maximum possible toxicity resistance score for ten different hyperpaths is plotted. To group similar results together, the score was rounded to the first digit. For some tradeoffs, all ten hyperpaths lead to the same or very similar toxicity resistance scores. For others (e.g. tradeoff 1 or 20), the results differ significantly.

Figure 2 .

 2 Figure 2.10: Maximum toxicity resistance score for hyperpaths where the score was omitted in the optimization problem. For each of the 47 tradeoffs, the maximum possible toxicity resistance score for ten different hyperpaths is plotted. These hyperpaths were computed without fixing the toxicity resistance score in the optimization problem. To group similar results together, the score was rounded to the first digit. Comparing the results to Figure 2.9, globally the scores are lower. Only two hyperpaths (tradeoff 19) have maximum scores that are over 20.

Figure 2 . 11 :

 211 Figure 2.11: Results of different hyperpaths for tradeoff 1. One thousand different hyperpaths were computed taking the toxicity resistance score into account (red points) and 1000 different hyperpaths were computed omitting it (blue points).

Figure 2 .

 2 Figure 2.12: Results of different hyperpaths for tradeoff 3.

Figure 2 .

 2 Figure 2.13: Results of different hyperpaths for tradeoff 6.

Figure 2 .

 2 Figure 2.14: Results of different hyperpaths for tradeoff 13.

Figure 2 .

 2 Figure 2.15: Results of different hyperpaths for tradeoff 47.

Figure 3 . 1 :

 31 Figure 3.1:Expected active reactions for different pulse experiments. These essential reactions along with their expected directions are highlighted in red whereas other non-essential reactions (but which nonetheless could be chosen) are depicted in grey. Each pulse is indicated by the short yellow arrow (Glc: glucose; Pyr: pyruvate and Suc: Succinate). During the glucose pulse, the glycolysis reactions (depicted in green) should be active in order to generate ATP from the hydrolysis of glucose. On the other hand, the pyruvate and succinate pulse experiments should show gluconeogenesis activation (also depicted in green but in the opposite sense), generating glucose-6-phosphate from these two carbon sources. Furthermore, the TCA cycle (depicted in blue) can be fed from pyruvate during the pyruvate and glucose pulses. During the succinate pulse, the overflow in the TCA cycle should lead to the production of pyruvate with a subsequent activation of gluconeogenesis to produce biomass precursors. The pentose phosphate pathway (depicted in purple) is most likely active in all pulses in order to generate biomass precursors; however, since this pathway is a mere interconversion of carbohydrates, there is no particular expectation as to the actual direction of these reactions.

Figure 3 . 2 :

 32 Figure 3.2: E. coli core model -Results for pyruvate pulse (λ = 0.1, = 5). The metabolites that were given as input are highlighted in blue if the corresponding input deltas were below zero and red if they were above zero. Reactions that are highlighted in orange were chosen in almost all of the enumerated solutions. Reactions that are yellow were chosen only in around half of the solutions. White reactions were not chosen in any solution. The expected reactions of the gluconeogenesis and part of the TCA cycle are active in all 100 solutions. The reversible reactions of the gluconeogenesis are chosen in the correct direction. The figure was created using Escher(King et al., 2015a).

  preferable. Besides the biomass precursor G3P, Totoro proposed the production of alpha-D-ribose 5-phosphate (R5P) via ribose-5-phosphate isomerase (RPI) and the production of D-erythrose 4phosphate (E4P) via Transketolase (TKT2). The results for λ = 0.1 and = 5 are shown in Figure3.2 (see Supplementary Figure A.5 for = 10).

. 3 .Figure 3 . 3 :

 333 Figure 3.3: E. coli core model -Results for glucose pulse.The labels for reaction and metabolite names were omitted. The metabolites that were given as input are highlighted in blue if the corresponding input deltas were below zero and red if they were above zero. Reactions that are highlighted in orange were chosen in almost all of the enumerated solutions. Reactions that are yellow were chosen only in around half of the solutions. White reactions were not chosen in any solution. In (a) and (b), even for λ = 0.1 the active reactions remain disconnected. In (c) and (d), lowering allowed for a fine-tuning of the solution which made it possible to obtain complete pathways. All figures were created using Escher(King et al., 2015a).

Figure 3 . 4 :

 34 Figure 3.4: E. coli core model -Results for succinate pulse (λ = 0.1, = 5). The metabolites that were given as input are highlighted in blue if the corresponding input deltas were below zero and red if they were above zero. Reactions that are highlighted in orange were chosen in almost all of the enumerated solutions. Reactions that are yellow were chosen only in around half of the solutions. White reactions were not chosen in any solution. The reactions of the gluconeogenesis pathway and the reactions that transform succinate in the TCA cycle and subsequently into pyruvate are active in all 100 solutions. The figure was created using Escher(King et al., 2015a).

Figure 3 . 5 :

 35 Figure 3.5: Carbon sources and closely connected amino acids. The production cost for individual amino acids can depend on the available carbon source, and reactions close to the entry point of the carbon source might have considerably higher fluxes. This figure shows the entry points for different carbon sources and closely connected amino acids. For example, Pyruvate is a direct precursor for valine production.

5 :

 5 Number of reaction occurrences in 19 different knockout sets.

Figure A. 2 :

 2 Figure A.2: E. coli core model -Results for pyruvate pulse (λ = 0.9, = 5). Similar reactions are active as for λ = 0.9 and = 10. An important difference is that the biomass reaction is part of the solution. The figure was created using Escher (King et al., 2015a).

Figure A. 4 :

 4 Figure A.4: E. coli core model -Results for pyruvate pulse (λ = 0.5, = 5). The results are similar to λ = 0.5 and = 10. However, again the biomass reaction is active in all solutions. The figure was created using Escher (King et al., 2015a).

Figure A. 6 :

 6 Figure A.6: E. coli core model -Results for glucose pulse (λ = 0.9, = 5). Similar to the results of the pyruvate pulse for λ = 0.9, the active reactions in the solutions are disconnected. The figure was created using Escher (King et al., 2015a).

Figure A. 10 :

 10 Figure A.10: E. coli core model -Results for glucose pulse (λ = 0.1, = 1.2).The previous result could be improved even further by decreasing λ to 1.2. Reactions of the glycolysis and the TCA cycle are active in all 100 solutions. The active reactions of the glycolysis pathway that are reversible are chosen in the correct direction. Furthermore, in contrast to the previous results, the biomass reaction is part of the solution. However, the reaction that transports glucose is only active in 64 out of 100 solutions. Lowering λ to 1.1 rendered the optimization problem infeasible. The figure was created using Escher(King et al., 2015a).

Figure A. 12 :

 12 Figure A.12: E. coli core model -Results for succinate pulse (λ = 0.9, = 10). We can see similarities to the results of the other pulse experiments. For λ = 0.9 the active reactions are disconnected. The figure was created using Escher(King et al., 2015a).

Figure A. 14 :

 14 Figure A.14: E. coli core model -Results for succinate pulse (λ = 0.5, = 10). The results are already more connected than for λ = 0.5 but we are not able to obtain the expected gluconeogenesis pathway. The figure was created using Escher (King et al., 2015a).

Figure A. 16 :

 16 Figure A.16: E. coli core model -Results for succinate pulse (λ = 0.1, = 10). The reactions of the gluconeogenesis pathway and the reactions that transform succinate in the TCA cycle and subsequently into pyruvate are active in all 100 solutions. The figure was created using Escher (King et al., 2015a).

  Metabolic networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Metabolic network reconstruction and databases . . . . . . . . . . . . . 1.2.2 Graph representation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Stoichiometric matrix . . . . . . . . . . . . . . . . . . . . . . . . . .

Contents 1.1 Linear programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Multi-objective linear programming . . . . . . . . . . . . . . . . . . . 1.2 1.2.4 Constrained-based modeling for metabolic networks . . . . . . . . . . . 1.3 Predicting optimal intervention strategies . . . . . . . . . . . . . . . . . . . . 1.4 Metabolic shifts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Theoretical background

  ). The metabolites S 1 , S 2 , S 3 and S 4 are the sources of the metabolic network. Metabolite T is the target node whose topological sources should be identified. For example, if M 0 = {S 1 , S 2 , S 3 }, it is possible to produce A because both substrates of reaction r 1 are present. Therefore, M 1 = {S 1 , S 2 , S 3 , A}. Now, both substrates of reaction r 2 are available and C can be produced which means that subsequently also T can be produced. Afterwards, no more changes can be made and a fixed point is reached. Hence, F wd({S 1 , S 2 , S 3

Table A .

 A

	Id Biomass Ethanol Toxicity
	1 0.131931 17.6169 37.0072
	2 0.138936 17.6120 36.0071
	3 0.145084 17.5654 35.0066
	4 0.157094 17.6384 26.0035
	6 0.165196 17.5216 26.0035
	9 0.180441 16.6031 37.0120
	11 0.187788 16.6035 36.0115
	13 0.185612 15.0541 43.0187
	19 0.103374 13.4770 52.0291
	47 0.347631 10.2221 43.0447

2. A reduced version that contains only the extreme points that will be referred to in the text is shown in 2.1.

Table 2 .

 2

1:

Computed tradeoffs between biomass production, ethanol production and toxicity score -short version. Only the extreme points in the Pareto front were computed. Lower bounds for biomass production was set to 0.1, lower bounds for ethanol production to 10. The list is sorted by descending ethanol production. The list presented here contains only those tradeoffs that are explicitly mentioned in the text. For a full list, we refer the reader to the Supplementary Table

A.2.

Table 2 .

 2

	Id Max. biomass Min. ethanol Max. ethanol Max. toxicity score
	1	0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.370	13.003
		0.1642	14.677	15.021	13.003
		0.1642	14.677	15.021	13.003
		0.1643	17.101	17.297	21.008
	2	0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1642	14.689	15.033	13.003
		0.1643	14.788	15.137	13.003
		0.1642	14.676	15.020	13.003
	4	0.1716	17.433	17.473	11.002
		0.1716	17.433	17.473	11.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
		0.1777	17.348	17.388	10.002
	9	0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
		0.2870	12.716	12.732	20.041
	11	0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041
		0.2874	12.717	12.733	22.041

2: Results for ten different hyperpaths for certain tradeoffs. For tradeoffs 1, 2, 4, 9 and 11 the production values are shown for ten different hyperpaths each.

Table 2 . 3 :

 23 H-Id Max. biomass Min. ethanol Max. toxicity score List of hyperpaths that were chosen for the random exploration. The first column indicates the id of the hyperpath. A total of 13 hyperpaths of tradeoff 1 were chosen based on their minimum ethanol productions and their maximum toxicity scores. The table is sorted by decreasing toxicity scores. The last hyperpaths (#10) has a slightly lower maximum toxicity score.

		799	0.1690	16.713	22.016
		800	0.1686	16.747	22.015
		656	0.1706	16.753	22.014
		657	0.1706	16.753	22.014
		599	0.1661	16.815	22.011
		31	0.1659	16.831	22.011
		537	0.1659	16.831	22.011
		536	0.1659	16.831	22.011
		638	0.1653	16.899	22.011
		169	0.1671	16.520	22.011
		879	0.1657	16.961	22.010
		638	0.1652	16.992	22.009
		10	0.1643	17.101	21.008
	10	1	0.1691	14.829	17.327	21.024
	10	1	0.1690	16.119	17.331	21.016
	10	2	0.1692	15.912	16.687	21.012
	169	1	0.1721	14.773	16.819	22.024
	536	1	0.1709	15.487	16.698	21.019
	536	2	0.1709	14.627	15.864	21.012
	536	2	0.1707	14.836	15.978	21.011
	537	1	0.1709	14.220	15.864	21.012
	537	1	0.1709	14.220	15.864	21.012
	537	1	0.1709	15.427	15.864	21.012
	537	1	0.1709	15.427	15.864	21.012
	599	1	0.1720	15.177	15.689	21.011
	638	1	0.1702	16.943	17.109	21.010
	638	2	0.1739	14.637	16.524	22.028
	638	2	0.1739	15.394	16.563	22.024
	656	1	0.1709	15.801	16.702	21.014
	658	2	0.1729	16.211	16.817	22.021
	879	1	0.1706	15.177	16.421	21.018

H-Id Run Size Max. biomass Min. ethanol Max. ethanol Max. toxicity score

Table 2 .

 2

	Gene associations
	Sbml Id Reaction name

4: Size and production values for selected knockout sets. In total, 18 interesting knockout sets got selected based on their size and their values for ethanol productions and toxicity scores. The smallest knockout sets that were enumerated contained 17 reactions.

Table 2 . 5 :

 25 Number of reaction occurrences in 19 different knockout sets -short version. The table shows how often a reaction occurred in the 19 selected knockout sets. The number of occurrences is shown in the last column. In this version, only reactions that appeared at least five times are shown. For the full version, see Supplementary Table A.5. In total, 74 reactions appeared in the 19 knockout sets. Problematically, quite a few of the proposed reactions are transports that might be difficult to knock out in vivo.

Table 3 . 1 :

 31 ). Comparison of different objective values for the best runs for each experiment. Since we are not fixing the objective value of the first solution in our optimization problem, the objective values for the subsequent solutions can be worse. In this table, we are comparing the difference in the objective values between the first solution and the 100th solution. In addition to the absolute differences, also the percentage of how much the objective value worsened compared to the first solution is displayed. The underlying optimization problem is a minimization problem. Therefore, smaller objective values are better.

	Pulse experiment	1st sol. 100th sol. Abs. diff.	%
	Pyruvate (λ = 0.1, = 5)	-32.1394	-30.6635	1.48	5.5%
	Glucose (λ = 0.1, = 1.2)	5.3830	6.5582	1.18 21.8%
	Succinate (λ = 0.1, = 5) -158.1770 -157.5760	0.60	0.4%

Table 3 . 2 :

 32 Results for pyruvate (λ = 0.1, = 1.0). In total, 10 solutions were computed. The table shows the objective value for each solution and how close this value is to the optimum (in %).

	Iteration Objective value Optimality (%) Runtime (h)
	1	-135.112	5.66
	2	-134.915	5.48
	3	-134.879	5.00
	4	-134.902	6.19
	5	-134.976	5.00
	6	-135.309	5.00
	7	-135.303	6.62
	8	-135.142	6.31
	9	-135.117	5.65
	10	-135.146	6.45

Table 3 . 3 :

 33

	1	-147.134	8.39
	2	-147.153	8.35
	3	-147.744	6.58
	4	-147.75	8.09
	5	-147.744	7.72
	6	-147.772	7.47
	7	-147.809	8.32
	8	-147.998	7.30
	9	-148.259	7.47
	10	-148.331	8.08

Theoretical background model reversible reactions, bidirectional hyperarcs can be used.

Conclusions and perspectivesHence, it might be interesting to see if our method is adaptable to also integrate, for instance, transcriptomic and/or proteomic data.

Acknowledgments

I would like to thank my advisors Marie-France Sagot and Susana Vinga and all other members of the Erable team; in particular Arnaud Mary for our discussions and my co-authors Ricardo Andrade and Mariana Ferrarini; special thanks to my colleagues Carol Moraga and Marianne Borderes who always provided support and comfort; my family for their support and for not asking too many questions.

List of Tables

The interval for each metabolite was calculated based on the baseline concentration and the fold change for the pseudo-steady state taking into account the given small variations for the baseline measurements.

To add glucose as a source, ∆ min glucose was set to -50. Oxygen was added as a source in the same way. Biomass was added as a sink. The interval for each metabolite was calculated based on the baseline concentration and the fold change for the pseudo-steady state taking into account the given small variations for the baseline measurements. To add glucose as a source, ∆ min succinate was set to -100. Oxygen was added as a source in the same way. Biomass was added as a sink. Higher values for sources and sinks were chosen due to larger calculated ∆ min and ∆ max than in the other two experiments.

Network model

Active reactions in the E. coli core model

The following figures show the results for the three pulse experiments for different parameters. The metabolites that were given as input are highlighted in blue if the corresponding input deltas were below zero and red if they were above zero. Reactions that are highlighted in orange were chosen in almost all of the enumerated solutions. Reactions that are yellow were chosen only in around half of the solutions. White reactions were not chosen in any solution.

Appendix

Figure A.1: E. coli core model -Results for pyruvate pulse (λ = 0.9, = 10). The active reactions are disconnected. Since λ = 0.9, the optimization prioritizes the minimization of the number of active reactions, fewer active reactions are thus chosen in total and the accumulation/depletion of unmeasured metabolites is higher. The figure was created using Escher (King et al., 2015a).