Emilie Allart

Joachim Niehren

Cristian Versari

Computing Difference Abstractions of Metabolic Networks Under Kinetic Constraints

Keywords: Gene knockout prediction, reaction networks, constraints, systems biology, synthetic biology, metabolism

Algorithms based on abstract interpretation were proposed recently for predicting changes of reaction networks with partial kinetic information. Their prediction precision, however, depends heavily on which heuristics are applied in order to add linear consequences of the steady state equations of the metabolic network. In this paper we ask the question whether such heuristics can be avoided while obtaining the highest possible precision. This leads us to the first algorithm for computing the difference abstractions of a linear equation system exactly without any approximation. This algorithm relies on the usage of elementary flux modes in a nontrivial manner, first-order definitions of the abstractions, and finite domain constraint solving.

Introduction

Flux balance analysis [START_REF] Orth | What is flux balance analysis?[END_REF][START_REF] Papin | Comparison of network-based pathway analysis methods[END_REF] can be used to predict the effect of influx changes of metabolic networks at steady state. Such predictions can be based on reasoning with linear equations systems that describe the rates of the reactions in a steady state of the metabolic network, by using Gaussian elimination, elementary flux modes (EFMs) [START_REF] Lotz | Elementary flux modes, flux balance analysis, and their application to plant metabolism[END_REF], or optimisation methods [START_REF] Maranas | Flux Balance Analysis and LP Problems[END_REF][START_REF] Facchetti | Partial inhibition and bilevel optimization in flux balance analysis[END_REF]. Most importantly, precise quantitative kinetic information is not required in contrast to classical mathematical analysis methods for reaction networks [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors[END_REF][START_REF] Calzone | BIOCHAM: an environment for modeling biological systems and formalizing experimental knowledge[END_REF][START_REF] Hoops | Copasi-a complex pathway simulator[END_REF][START_REF] Fages | Inferring reaction systems from ordinary differential equations[END_REF]. In fact, even when the kinetic functions associated to chemical reactions are known, the values of rate constants are most often missing, since it is difficult to measure them experimentally in the precise state of the regulation of the metabolic network at the time point of interest.

Recently, abstract interpretation [START_REF] Cousot | Systematic design of program analysis frameworks[END_REF][START_REF] Fages | Abstract interpretation and types for systems biology[END_REF][START_REF] Danos | Abstracting the differential semantics of rule-based models: Exact and automated model reduction[END_REF] has been exploited to design novel algorithms [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF][START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF][START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF]] that can use partial kinetic information beneficially for predicting changes of metabolic networks. They can in particular exploit the knowledge about the enzymes and inhibitors. Similarly to flux balance analysis, the linear equations describing steady states are used, but in addition to them, kinetic constraints are inferred from the partial kinetic information of inhibitors and enzymes.

The steady state equations and the kinetic constraints enable gene knockout predictions based on abstract interpretation (in the finite relational structure ∆ 6 ), based on the linear equations from the metabolic network and the constraints on its regulatory control. The unknown kinetic parameters are abstracted away, by interpretation over some finite relational structure, that contains a finite number of abstract differences rather than concrete differences in R 2 + . Eventually, the prediction algorithm will apply a finite domain constraint solver for ∆ 6 that we implemented in Minizinc [START_REF] Rendl | Minisearch: A solver-independent meta-search language for minizinc[END_REF] to enumerate all the changes that may or must lead to the target change.

The prediction quality of abstract interpretation approaches heavily depends on heuristics that find and add linear equations entailed by the steady state equations before constraint solving over ∆ 6 . This is necessary to enable global reasoning, since local reasoning alone is not able to Fig. 1. A toy metabolic network with a simple cycle. Fig. 2. A glimpse of the formal model from [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] of leucine production in B. subtilis.

deal with cycles in metabolic part of the network as we will illustrate in Section 2. On the other hand, it is impossible to add the infinity of all entailed linear equations before abstract interpretation. Therefore, these algorithms can at best approximate the abstraction of differences of solution set of linear equations. Whether this abstraction can be computed exactly is a long standing open problem, as well as how to measure the quality of approximation heuristics.

In this paper, we present the first exact algorithm that can compute the ∆ 6 difference abstraction of the solution set of a linear equation system without any overapproximation. We apply it to the prediction of leucine overproduction, a benchmark task that is best studied with abstract interpretation. In this case, we need to deal with kinetic constraints in addition that are naturally interpreted over ∆ 6 . It turns out that a new heuristic that we also propose in the present paper does indeed compute well the difference abstraction at this benchmark task from systems biology, although being inexact in the general case. The main advantage of this heuristic is that it outperforms the exact algorithm dramatically in computation time: only 5 minutes are needed for the knockout prediction rather than 5 hours with the always exact algorithm.

Qualitative reasoning on metabolic reaction networks

The application of abstract interpretation to the analysis of metabolic reaction networks is based on an intuitive qualitative reasoning. Its aim is to predict how a living organism or its environment should be changed in order to maximize the production of some metabolite of interest, without exact knowledge of the quantitative parameters of the system. A change of the organism is represented for example by any modification of its genome, such as a gene knockout. A change of the environment is typically represented by the modification of the culture medium, which results in an increase or decrease of some inflows. Since changes can be arbitrarily combined to obtain or improve the wanted results, the problem that we tackle is highly combinatorial.

Examples of formal metabolic networks on which abstract interpretation can be applied are shown in Fig. 1 and Fig. 2. Fig. 1 shows a toy metabolic network with a simple cycle that will be used in the following to introduce the key ideas of the reasoning. In this network, metabolites are displayed in yellow rounded boxes, while reactions are in gray squared boxes, in the tradition of Petri nets. Reactions with dotted contour are inflows or outflows of the system. Fig. 2 shows a glimpse of a bigger metabolic reaction network with regulation, where the regulatory part is represented by the enzymes in blue rounded boxes. The full network of Fig. 2 models a part of the metabolism of the gram positive bacterium B. subtilis. Our benchmark application -taken from [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] -is the overproduction of one of the metabolites of this network, the branched chain amino acid Leucine (Leu). This amino acid is a precursor of surfactin, a non ribosomal peptide with several applications in food and pharmaceutical industry.

Let us now reconsider the toy model with a simple cycle in Fig. 1. This network is composed of two chemical species A and B. The species A is continuously produced by an inflow at a fixed rate X, and is transformed into B by the reaction with rate V . The inverse reaction with rate V -transforms B back into A. The species B has an outflow with rate Y . All reactions but the inflow are controlled internally by the system. The outflow rate Y in particular is determined by the concentration of B which in turn depends on the rates V and V -.

The only possible change in this toy network is the increase or decrease of the inflow X. In order to illustrate the reasoning method, we set the increase of Y as our final target. As usually done in flux balance analysis [START_REF] Orth | What is flux balance analysis?[END_REF][START_REF] Papin | Comparison of network-based pathway analysis methods[END_REF], we consider the system at the steady state, that gives us the following linear system of equations:

∃V ∃V -. V = X + V -∧ V = Y + V - (1)
The existential quantifiers for V and V -allow us to hide the internal behaviour of the network, so to project to inflows and outflows. While the consequences of steady state equations may be difficult to interpret without using Gauss' algorithm, in this particular case it is easy to see that by subtracting the first equation in (1) from the second we obtain an important relation:

X = Y (2)
Eq. ( 2) tells us that the only way for Y to increase is that X increases too. We formalize now this intuition by means of abstract interpretation, this time applied to concrete differences in R 2 + . Concrete differences capture the essence of a change in the spirit of [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF]: a change of the value of X for instance can be thought of as a pair of positive reals (r before , r after ) representing the value of X at steady state respectively before and after the modification of the environment. We need to consider positive reals since the rates of irreversible reactions are positive. We call the above concrete differences an increase if r before < r after , a decrease if r before > r after and a nochange if r before = r after . This intuition motivates the usage of abstract values in ∆ 3 = { , , ∼ ∼ ∼ } where = (0, 1) represents an increase, = (1, 0) a decrease, and ∼ ∼ ∼ = (0, 0) a no-change. The canonical mapping of concrete differences in R 2 + to the abstract differences in ∆ 3 can be seen as a homomorphism between the relational structures R 2 + to ∆ 3 . This algebraic view of abstractions as homomorphisms enables various generalizations. An example is the abstraction from R 2 + to ∆ 6 -as considered for gene knockout prediction [START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF][START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] -which refines each of the three R 2 + equivalence classes produced by ∆ 3 into two parts in ∆ 6 .

Any abstraction of concrete differences enables some form of abstract qualitative reasoning [START_REF] Forbus | Qualitative reasoning[END_REF] based on operations of the relational structure of abstract differences, that can be used for change prediction in systems biology. Let us illustrate how to reason with ∆ 3 . As a first example, assume that we know for some reason that X and V -both increase, that is X = and V -= . Then we can use the first equation in [START_REF] Allart | Computing sign abstractions of linear systems[END_REF] to deduce for sure V will also increase, since + ∆3 = . The full table defining the summation operator + ∆3 on abstract difference is given in Fig. 3.

The above qualitative reasoning method, however, is quite weak when relying only on the steady state equations in [START_REF] Allart | Computing sign abstractions of linear systems[END_REF]. The main reason is that all reasoning steps are local, so that they overlook global properties of the network that are arising for example with metabolic cycles. To see this, suppose that we want to predict which environmental change may lead to an increase of Y . We can use the second equation in [START_REF] Allart | Computing sign abstractions of linear systems[END_REF] to infer some constraints on the values of V and V -: if Y = , we can for example infer that V cannot be if V -is . In fact, + can never be equal to . However, V = , V -= is a partial solution that seems compatible with an increase of Y . This illustrates that we cannot infer any constraint on X with this kind of local reasoning.

Equation X = Y in Eq. ( 2), in contrast, expresses a global property of the network that immediately implies that the only value for X compatible with an increase of Y is . In other words, when reasoning with equations over abstract differences, Eq. ( 2) is no longer an implicit consequence of the system [START_REF] Allart | Computing sign abstractions of linear systems[END_REF]. Therefore, to be taken into account, such equations must be explicitly included in the system before applying the abstract reasoning. Unfortunately, the number of entailed linear equations is infinite in general. For instance, in our small network all the equations of the family nV + mX = nV -+ (m + 1)Y for any two naturals n, m are consequences of (1). So, instead of trying to infer the set of all the consequences of our system, we may try to compute a "good" subset of it, by including only the consequences that more heavily constrain the variables. The prediction quality of the existing approaches heaviliy depends on the heuristics chosen for adding entailed linear consequences.

One evident advantage of Eq. ( 2) is its small number of variables (there are only two: X and Y ). If we consider an individual linear equation, the intuition is that removing a variable increases the constraining power on the remaining variables. So, we propose as a heuristics the inclusion of all linear consequences involving a minimal subset of variables. This idea is at the core of the first result that we present in this paper: a heuristic algorithm that enriches the steady state equations of the metabolic network with minimal support consequences before applying abstract interpretation. While this algorithm has been internally used for some time to increase the precision of the abstract interpretation, some key questions have been always open about it:

1. is the set of minimal support equations complete? that is, does it represent all the deducible constraints on the abstract system? 2. are these constraints sufficient to compute the exact set of abstract solutions? 3. if this is not the case, is there a method to compute them?

The definitive answers to these questions is the second main contribution of this paper. We show in particular that the above heuristic does not cover all the entailed constraints on the abstract system, i.e. it does not allow to compute the exact abstraction of a linear system in general. Intuitively, this happens because the approach takes for granted that the abstract reasoning is based on the linear system computed at the steady state, that is on the matrix equation AX = 0 where A is the stoichiometry matrix associated to the metabolic network, and X is the set of metabolic flows representing the unknowns of the linear system. However, it is easy to notice that as soon as concrete differences are introduced in the reasoning, there is not only one, but actually two linear systems to consider: one before the environmental change, and one after it, that is one system for each value of the pairs representing concrete changes. Informally, we should therefore consider a bigger matrix equation including somehow both AX before = 0 and AX after = 0.

This idea is the starting point of development of the main contribution of the present paper: a method for the exact computation of the abstraction of a linear system, that we call the exact algorithm. This method not only provides us with the counterexamples to the exactness of the heuristic based on minimal support consequences, but gives us also an exact measure of its goodness as well as of the goodness of all the other heuristics used to improve our abstract analysis. Remarkably, both the heuristic and the exact algorithm have their root in the rewriting of a linear system in terms of its EFMs. The key difference between the two methods lies in the choice of the linear system (i.e. the matrix equation) initially used to compute the EFMs.

Preliminaries

Set Notation.We start with usual notation for sets. Let N be the set of natural numbers and R + the set of positive real numbers, both including 0. For any set A and n ∈ N, the set of n-tuples of elements in A is denoted by A n . The i-th projection function on n-tuples of elements in A, where 1 ≤ i ≤ n is the function π

i : A n → A such that π i (a 1 , . . . , a n ) = a i for all a 1 , . . . , a n ∈ A. If A is finite the number of elements of A is denote by |A|.
Σ-Algebras and Σ-Structures. We next recall the notions of Σ-algebras, Σ-structures, and homomorphism between Σ-structures. Let Σ = ∪ n≥0 F (n) C be a ranked signature. The elements of f ∈ F (n) are called the n-ary function symbols of Σ and the elements in c ∈ C its constants. Definition 1. A Σ-algebra S = (dom(S), . S ) consists of a set dom(S) and an interpretation . S such that c S ∈ dom(S) for all c ∈ C, and f S : dom(S) n → dom(S) for all f ∈ F (n) and n ∈ N.

We next reinterpret n-ary function symbols of Σ as n+1-ary relation symbols, so that we can reuse the same signature Σ for defining Σ-structures. Definition 2. A Σ-structure ∆ = (dom(∆), . ∆ ) consists of a set dom(∆) and an interpretation . ∆ such that c ∆ ∈ dom(∆) for all c ∈ C and f ∆ ⊆ dom(∆) n+1 for all f ∈ F (n) and n ∈ N.

In this manner, any Σ-algebra is also a Σ-structure since any n-ary function is an n + 1-ary relation. Note also that symbols in F (0) are interpreted as monadic relations in Σ-structures, i.e., as subsets of the domain, in contrast to constants in C that are interpreted as elements of the domain. Definition 3. A homomorphism between two Σ-structures S and ∆ is a function h : dom(S) → dom(∆) such that for c ∈ C, n ∈ N, f ∈ F (n) , and s 1 , . . . , s n+1 ∈ dom(S):

1. h(c S ) = c ∆ , and 2. if (s 1 , . . . , s n+1 ) ∈ f S then (h(s 1 ), . . . , h(s n+1 )) ∈ f ∆ .
If we consider n + 1-ary relations as n-ary set valued functions, the second condition can be rewritten equivalently as h(f

S (s 1 , . . . , s n )) ⊆ f ∆ (h(s 1 ), . . . , h(s n )). For Σ-algebras, this condi- tion is equivalent to h(f S (s 1 , . . . , s n )) = f ∆ (h(s 1 ), . . . , h(s n )).

Σ-Abstractions

Throughout the paper we will use the signature Σ = F (2) C with two binary function symbols in F (2) = {+, * } and two constants C = {0, 1}. We will only consider Σ-algebras S in which + S and * S are associative and commutative, with neutral element 0 S and 1 S respectively. Example 4. The set of positive real numbers R + can be turned into a Σ-algebra with domain R + , by interpreting + as the addition of positive real numbers + R+ , * as the multiplication of positive real numbers * R+ , and interpreting the constants by themselves 0 R+ = 0 and 1 R+ = 1. We will deliberatly confuse the set R + with the Σ-algebra (R + , . R+ ).

Example 5. The set of Booleans B = {0, 1} ⊆ R + can be turned into a Σ-algebra with domain B by interpreting + B = ∨ B as disjunction, * B = ∧ B as conjunction, and the constants by themselves 0 B = 0 and 1 B = 1. We will deliberatly confuse the set B with the Σ-algebra (B, . B ).

We can abstract positive real numbers into booleans by defining a function h

B : R + → B such that h B (0) = 0 and h B (r) = 1 for all r ∈ R + \ {0}. Lemma 6. The function h B : R + → B is a homomorphism between Σ-algebras.
The homomorphism h B is the prime example of what we will call a Σ-abstraction. Definition 7. A Σ-abstraction is a homomorphism between Σ-structures S and ∆ such that dom(∆) ⊆ dom(S).

Abstracting Concrete Differences

Concrete differences are pairs of positive real numbers such as (r before , r after ) ∈ R 2 + in the example Section 2. We show how to abstract concrete differences into abstract differences in some finite Σ-structure. The Tuple Σ-Algebra S n . For any Σ-algebra S and natural number n ∈ N we define the Σalgebra of n-tuples S n = (dom(S) n , . S n ) such that for all s 1 , . . . , s n , s 1 , . . . , s n ∈ dom(S) and ∈ F (2) : (s 1 , . . . , s n ) S n (s 1 , . . . , s n ) = (s 1 S s 1 , . . . , s n S s n )

The constants c ∈ C are interpreted as c S n = (c S , . . . , c S ). Note that if 0 S is the neutral element of + S , then 0 S n is the also the neutral element of + S n . In analogy, if 1 S is the neutral element of * S then 1 S n is also the neutral element of * S n . Furthermore, the associativity and commutativity of + S n and * S n inherit from + S and * S respectively. Note that we deliberatly confuse the set R 2 + with the Σ-algebra (R 2 + , . R 2 + ) with our notation. Given this, it follows from the above, that the algebra R 2 + has the neutral element (0, 0) for + R 2 + and the neutral element (1, 1) for * R 2 + , and that these operations are associative and commutative. For any function h : A → B and n ∈ N we define the function h n : A n → B n such that h n (a 1 , . . . , a n ) = (h(a 1 ), . . . , h(a n )) for all a 1 , . . . , a n ∈ A.

Lemma 8. If h is Σ-abstraction from S to ∆ then h n is a Σ-abstraction from S n to ∆ n .
Abstractions of Concrete Differences. A generic manner to abstract concrete differences in R 2 + is to start with a finite set ∆ ⊆ R 2 + of so called abstract differences, and some function h : R 2 + → ∆ that says how to abstract any concrete differences to some abstract difference. The function h defines a partition of R 2 + into the equivalences classes of concrete differences that are mapped to the same abstract difference.

Given such a function h, there is a unique manner to define an interpretation . ∆ such that (∆, . ∆ ) becomes Σ-structure and h a Σ-abstraction. For any constant c ∈ C we have to define c ∆ = h(c R 2 + ) and for any function symbol ∈ F (2) we have to define a ternary relation ∆ , which seen as set-valued function ∆ : ∆×∆ → 2 ∆ must satisfy for all abstract values δ 1 , δ 2 ∈ ∆:

δ 1 ∆ δ 2 = {h(r 1 R+ r 2 , r 1 R+ r 2 ) | h(r 1 , r 1 ) = δ 1 , h(r 2 , r 2 ) = δ 2 } Lemma 9. h : R 2 + → ∆ is a Σ-abstraction.
The Σ-Structure ∆ 3 . Our next objective is to recall the abstraction of concrete differences into the finite Σ-structure with domain ∆ 3 = { , , ∼ ∼ ∼ } that is well-known from qualitative reasonning (see e.g. [START_REF] Forbus | Qualitative reasoning[END_REF]). For this we start with the function h ∆3 (r, r ) ∈ ∆ 3 such that for any all r, r ∈ R + :

h ∆3 (r, r ) =    = (1, 0) if r > r = (0, 1) if r < r ∼ ∼ ∼ = (0, 0) if r = r
The relations + ∆3 and * ∆3 are is the symmetric closure of the relation in Fig. 3. Furthermore, h ∆3 : R 2 + → ∆ 3 is a Σ-abstraction by Lemma 9. The Σ-Structure ∆ 6 . We next recall the abstraction of concrete differences to the finite Σstructure with domain ∆ 6 = {↑, ↓, ∼, ⇑, ⇓, ≈} that was introduced for gene knockout prediction in [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF]. For defining this Σ-structure, we start with the function h ∆6 : R 2 + → ∆ 6 such that for any two numbers r, r ∈ R + :

h ∆6 (r, r ) =    ↑ = (1, 2) if 0 = r < r ↓ = (2, 1) if r > r = 0 ∼ = (1, 1) if r = r = 0 h ∆6 (r, r ) =    ⇑ = (0, 2) if 0 = r < r ⇓ = (2, 0) if r > r = 0 ≈ = (0, 0) if r = r = 0 δ δ δ + ∆ 3 δ δ * ∆ 3 δ { } { } { , ∼ ∼ ∼ , } { , ∼ ∼ ∼ , } ∼ ∼ ∼ { } { , ∼ ∼ ∼ } δ δ δ + ∆ 3 δ δ * ∆ 3 δ ∼ ∼ ∼ ∼ ∼ ∼ { ∼ ∼ ∼ } { ∼ ∼ ∼ } { } { } ∼ ∼ ∼ { } { , ∼ ∼ ∼ } c c ∆ 3 0 ∼ ∼ ∼ 1 ∼ ∼ ∼ Fig. 3. Interpretation of Σ-structure ∆ 3 . δ δ δ + ∆ 6 δ δ * ∆ 6 δ ↑ ↑ {↑} {↑} ↑ ↓ {↑, ∼, ↓} {↑, ∼, ↓} ↑ ∼ {↑} {↑} ↑ ⇑ {↑} {⇑} ↑ ⇓ {↑, ↓, ∼} {⇓} ↑ ≈ {↑} {≈} ⇑ ↓ {↑, ∼, ↓} {⇑} δ δ δ + ∆ 6 δ δ * ∆ 6 δ ⇑ ∼ {↑} {⇑} ⇑ ⇑ {⇑} {⇑} ⇑ ⇓ {↑, ∼, ↓} {≈} ⇑ ≈ {⇑} {≈} ∼ ∼ {∼} {∼} ∼ ≈ {∼} {≈} ∼ ↓ {↓} {↓} δ δ δ + ∆ 6 δ δ * ∆ 6 δ ∼ ⇓ {↓} {⇓} ≈ ≈ {≈} {≈} ≈ ↓ {↓} {⇓} ≈ ⇓ {⇓} {⇓} ↓ ↓ {↓} {↓} ↓ ⇓ {↓} {⇓} ⇓ ⇓ {⇓} {⇓} c c ∆ 6 0 ≈ 1 ∼ Fig. 4. Interpretation of Σ-structure ∆ 6 .
The relations + ∆6 and * ∆6 are the symmetric closure of the relations in Fig. 4. By Lemma 9, h ∆6 : R 2 + → ∆ 6 is a Σ-abstraction.

First-Order Logic

We first recall the standard first-order logic and then show how to enhance it with n-tuples without increasing the expressiveness. We fix a set of variables V (for instance V = N). The variables in V will be ranged over by x and y. The set of first-order expressions e ∈ E Σ and first-order formulas φ ∈ F Σ are constructed according to the abstract syntax in Fig. 5 from the symbols in the signature Σ, the variables in V, the first-order connectives, and the equality symbol . =. As shortcuts, we define the formula true = def 1 . =1 and for any sequence of formulas φ 1 , . . . , φ n we define ∧ n i=1 φ i as φ 1 ∧ . . . ∧ φ n which is equal to true if n = 0. We define formulas e . =0 by ¬e .

=0.

The semantics of a formula φ ∈ F Σ is a truth value, which depends on the Σ-structures S of interpretation and on a variable assignment α : V → dom(S). Any Σ-expressions e ∈ E Σ denotes a subset of values in dom(S), which will be singleton in case that S was a Σ-algebra. The semantic of equations e . =e is, as expected when interpreted over Σ-algebras S: the unique values of e and e in S must be equal. However, we will also need to interpret equations e . =e over Σ-structures. This is why, any expression e denotes a subset of the Σ-structure, not just a single element. We can then interpret equality as nondisjointness, i.e., e . =e holds in a Σ-structure S if e and e are interpreted as nondisjoint subsets of dom(S).

A variable assignment into a Σ-structure S is a partial function α : V → dom(S) for some subset V ⊆ V. Let S be a Σ-structure and α a variable assignment to S. Any Σ-expression e with V(e) ⊆ V can be interpreted as an element of dom(S) and any Σ-formula φ ∈ F Σ with V(φ) ⊆ V as a Boolean value. The set of solutions of a formula φ ∈ F Σ over a Σ-structure S with respect to some set of variables V ⊇ V(φ) is defined by:

sol S V (φ)={α : V → dom(S) | φ S,α = 1} If V = V(φ)
then we omit the index V , i.e., sol S (φ) = sol S V (φ). We next extend the first-order logic to n-tuples where the parameter n is fixed. In applications, we will use the case n = 2, that is the first-order logic with pairs. Back and forth compilers from first-order logic with and without tuples will be convenient later on.

The syntax of first-order logic with n-tuples is given in Fig. 6. The expressions o ∈ O n Σ are like the expression e ∈ E Σ except that variables x are now replaced by projection expressions

e . =e S,α = 1 if e S,α ∩ e S,α = ∅ 0 else φ ∧ φ S,α = φ S,α ∧ B φ S,α ¬φ S,α = ¬ B ( φ S,α ) ∃x.φ S,α = 1 if exists s ∈ dom(S). φ S,α[x/s] = 1 0 else
Fig. 5. Syntax and semantics of expressions and formulas of first-order logic.

o ∈ O n Σ ::= . π i(x) | c | o o where ∈ F (2) , c ∈ C, 1 ≤ i ≤ n ψ ∈ F n Σ ::= o . =o | ∃x.ψ | ψ ∧ ψ | ¬ψ where x ∈ V
Fig. 6. Σ-expressions and Σ-formulas of the first-order logic with n-tuples.

.

π i (x) where 1 ≤ i ≤ n.
The reason is that any variable does now denote an n-tuple of values, rather than a single value (while the interpretation of constants and function symbols remain unchanged). The only change in the semantics is that variables assignment β do now map to n-tuples of values of the domain, and that

. π i (x) S,β = {π i (β(x))}.
The set of solutions of a formula ψ ∈ F n Σ over a Σ-structure S is defined as follows:

n-sol S (ψ)={β : V(ψ) → dom(S) n | ψ S,β = 1}
We next show how to express any first-order formulas in F Σ , interpreted over a tuple algebra S n , by some formulas in F n Σ , interpreted over S. In a first step, we convert first-order expression in e ∈ E Σ -that we will interpret over the Σ-algebra S n -to n projected expressions Π i (e) ∈ O n Σ where 1 ≤ i ≤ n. For all operators ∈ F (2) and constants c ∈ C we define:

Π i (e e ) = def Π i (e) Π i (e ) Π i (x) = def . π i (x) Π i (c) = def c
In the second step, we convert any formula φ ∈ F Σ without tuples -that will be interpreted in the tuple algebra S n -to some formula φ n ∈ F n Σ with tuples.

e . =e n = def ∧ n i=1 Π i (e) . =Π i (e ) φ ∧ φ n = def φ n ∧ φ n ¬φ n = def ¬ φ n ∃x.φ n = def ∃x. φ n
Proposition 10. For any φ ∈ F Σ , Σ-structure S, and n ≥ 1:

sol S n (φ) = n-sol S ( φ n ).
Example 11. Let 3 = def 1 + 1 + 1 and 4 = def 1 + 1 + 1 + 1. The formula φ ∈ F Σ equal to:

3 * x + 4 * y . =0
then has the same solutions over R 2 + than the formula φ 2 ∈ F 2 Σ over R + below:

3 * . π 1 (x) + 4 * . π 1 (y) . =0 ∧ 3 * . π 2 (x) + 4 * . π 2 (

y)

.

=0

We next show how to rewrite any first-order formulas with tuples ψ ∈ F n Σ into some first-order formula ν(ψ) ∈ F Σ without tuples. The idea is to replace all projections π i (x) by new variables ν i (x). For this, we first fix n generators of fresh variables ν 1 , . . ., ν n : V → V. Second, we map any expression o ∈ O n Σ with projections to some expressions ν(o) ∈ E Σ without new variables:

ν( . π i (x)) = def ν i (x), ν(c) = def c, ν(o 1 o 2 ) = def ν(o 1 ) ν(o 2 ).
Third, we map any formula ψ ∈ F n Σ with projections to some formula ν(ψ) ∈ F Σ with fresh variables:

ν(o = o ) = def ν(o) = ν(o ) ν(¬ψ) = def ¬ν(ψ) ν(ψ ∧ ψ ) = def ν(ψ) ∧ ν(ψ ) ν(∃x.ψ) = def ∃ν 1 (x) . . . ∃ν n (x). ν(ψ)
Given an variable assignment β : V → dom(S) n with V ⊆ V, we define ν(β) : n i=1 ν i (V ) → dom(S) such that for all x ∈ V :

ν(β)(ν i (x)) = π i (β(x)))
Function ν is a bijection with range {α | α : n i=1 ν i (V ) → dom(S)}. The inverse of this function satisfies ν -1 (α)(x) = (α(ν 1 (x)), . . . , α(ν n (x)) for all α in the range and all x ∈ V . Proposition 12. For any ψ ∈ F n Σ , Σ-structure S, and n ≥ 1: n-sol S (ψ) = ν -1 (sol S (ν(ψ))).

We finish this section with a commutation property of the operator ν -1 with Σ-abstractions.

Proposition 13. For any subset R of variable assignments of type V → dom(S) where V ⊆ V, n ≥ 1, and Σ-abstraction h :

S → ∆: ν -1 (h • R) = h n • ν -1 (R).

Difference Abstraction

We next show how to recast the notions of difference abstractions from [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF][START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF][START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF] by applying our notion of Σ-abstractions to the Σ-algebra R 2 + . Let S be a Σ-algebra and V ⊆ V a subset of variables. For any two variable assignments α, α : V → dom(S), we define an assignment of variables to pairs of elements in the domain of the structure diff(α, α ) : V → dom(S) 2 -that we call the differences of α and α -such that for all variables x ∈ V , diff(α, α )(x) = (α(x), α (x)). For any subset R of variable assignments of type V → dom(S) we define the set of differences of assignments in R by:

diff(R) = {diff(α, α ) | α, α ∈ R}
Furthermore, for any Σ-abstraction h : S 2 → ∆ and subset R of difference abstractions of type V → dom(S) 2 we define the application of the abstraction h to R by h

• R = def {h • β | β ∈ R } Definition 14.
For any Σ-abstraction h : S 2 → ∆ and formula φ ∈ F Σ we define the difference abstraction of the S-solution set of φ by: sol S (φ

) ∆ = h • diff(sol S (φ))).
The original definition of sol(φ) ∆6 in [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF] did not make explicit the roles of diff and h ∆6 : R 2 + → ∆ 6 . Having done so, we can now see that the difference abstraction of the R + -solution sets of a formula is the R 2 + -solution set of the same formula.

Lemma 15. For any formula φ ∈ F Σ and Σ-structure S: diff(sol S (φ)) = sol S 2 (φ).

As an immediate consequence, we have for any Σ-abstraction h : S 2 → ∆ that sol(φ) ∆ = h ∆ • sol S 2 (φ). Our next objective is to show that we can overapproximate the set sol(φ) ∆ by sol ∆ (φ) (Corollary 19). In order to show this, let h :S → ∆ be a Σ-abstraction and α be a variable assignment into dom(S ): Lemma 16. For any expression e ∈ E Σ with V (e) ⊆ dom(α): h ( e S ,α ) ⊆ e ∆,h •α . Proposition 17. For any positive formula φ ∈ F Σ with V (φ) ⊆ dom(α): φ S ,α ≤ φ ∆,h •α . Theorem 18. For any positive formula φ ∈ F Σ : h • sol S (φ) ⊆ sol ∆ (φ).

Corollary 19. For any Σ-abstraction h : S 2 → ∆ and positive first-order formula φ ∈ F Σ :

sol S (φ) ∆ = h • diff(sol S (φ)) ⊆ sol ∆ (φ)
This is an obvious consequence from Theorem 18 and Proposition 15. If ∆ is finite then the set sol ∆ (φ) is finite. If furthermore φ is a conjunctive formula, we can therefore compute the set sol ∆ (φ) by a finite domain constraint solver (such as e.g. Minizinc [START_REF] Rendl | Minisearch: A solver-independent meta-search language for minizinc[END_REF]). In contrast, it remains unclear how to compute the finite set h • diff(sol S (φ)) for infinite structures S. The problem is open, even if φ is a system of homogenous linear equations and S = R + , so that the infinite set sol S (φ) has a finite solved form by a triangular matrix. This is the core of the objective that we tackle in the remainder of the present paper.

Objective

We formalize the full algorithmic problem that we will solve in this paper and illustrate its relevance to our benchmark application to systems biology.

Once having fixed the parameter ∆ ∈ {∆ 3 , ∆ 6 } the algorithmic problem has three inputs :

Linear system over R + : a first-order formula φ ∈ F Σ that represents a linear equation system.

(This formula typically captures the steady state equations of the model.) Constraint over ∆: a first-order formula φ ∈ F Σ∪∆ where the signature Σ is extended with additional constants of ∆ that will be interpreted by themselves. (This formula typically expresses the partial kinetic knowledge on the reactions in the model and the change target of the prediction task (e.g. overproduction of some metabolites).) Set of observable variables: a finite subset of variables V ⊆ V(φ) ∪ V(φ ). (This set typically contains the control variables such as inflows and gene knockouts as well as the target variables, but not the variables for the rate of the internal metabolic reactions. Since the number of solutions may be of cardinality |∆| |V | , it is essential to choose V as small as possible.)

The algorithmic output that has to be produced is the ∆-abstraction of differences of R + -solutions of φ, but constrained to the solutions of φ over the structure ∆, and projected to the variables of V . In other words, the algorithm will compute the following finite set where V = V(φ) ∪ V(φ ):

{β |V | β ∈ sol R+ V (φ) ∆ ∩ sol ∆ V (φ )}
The only restriction on the inputs that we will impose is that the first formula φ must represent a homogeneous system of linear equations in F Σ . For instance, the linear equation x -2y = 0 is captured by the equation x .

=y + y in F Σ where we cannot use the minus operator. See Section 9 for the general definition. The constraint φ ∈ F Σ∪∆ in contrast may be arbitrary, including nonlinear equations and universal quantifiers but must be interpreted abstractly over ∆, while the linear equation system is valid over R + . Note however, that any universal quantifiers in φ can be expressed by a simple conjunction, given that the interpretation domain ∆ is finite.

In Fig. 7 we illustrate how the inputs will be instantiated for our benchmark application of leucine overproduction (a glimpse of the reaction network was given in Fig. 2). In this case, we choose the parameter ∆ = ∆ 6 . The observable variables in stand for the rates of the inflows Threonine (x Thr ), Akb(x Akb ), etc, the rate of the target outflow Leucine (y Leu ), and the possible gene knockouts. The system of linear equations φ contains the steady state equations for the metabolic reactions in the network. These require that all metabolites must be produced and consumed at the same rate. For instance, Pyruvate is produced by the inflow of Threonine at rate x Thr and consumed by reactions 27 and 30 at rates r 27 and r 30 respectively. The yields the linear equation (Thr) of Fig. 7. Species Akb is produced by reaction 41 and consumed by 27, leading to the steady state equation (Akb). Leucine is produced by reaction 45 and consumed by its outflow, leading to equation (Leu). The constraint φ ∈ F Σ∪∆6 contains the overproduction target y Leu = ↑ in (target) and the kinetic constraints for all reactions, of which we show only constraint (27) for reaction 27. The kinetic constraints must be interpreted abstractly over ∆ 6 according the formal semantics of the modeling language [START_REF] Niehren | Predicting Changes of Reaction Networks with Partial Kinetic Information[END_REF] rather than concretely over R + . Therefore, the meaning of the constraints is purely qualitative and not at all quantitative. For instance, the constraint (27) states (beside others) that rate of reaction 27 increases if either of the concentrations of the reactants Pyr and Akb or of the enzymes IlvD, IlvBH or IlvC increase. Nothing is said about quantities of these increases.

Exact Algorithms

We now present an exact solution of the problem presented in the previous section. Our approach is to characterize the abstraction of the solution set of a linear equation system as the solution set of some first-order formula over the abstract domain. We consider the abstractions h B , h ∆3 , and h ∆6 in this order.

Characterizing B-Abstractions. We now present a result from [START_REF] Allart | Computing sign abstractions of linear systems[END_REF] that shows that the boolean abstraction of the R + -solution set of a mixed linear system can be computed exactly. The development of this result was motivated by the needs of the present paper, but given that it is of independent interest and nontrivial to prove, we decided to present it independently.

Any natural numbers n can be described by the expression n = def n i=1 1 in E Σ . This permits to define linear equations as equations in F Σ that have the form:

n 1 * x 1 + . . . n m * x m . =n m+1 * x m+1 + . . . + n p * x p (3)
where m, p, n 1 , . . . , n p ∈ N and x 1 , . . . , x p ∈ V. Mixed Linear Systems. A product-zero-equation in F Σ is an equation of the form x * y . =0 where x, y ∈ V. A mixed linear system is a conjunctive formula in F Σ of the form ∃z. L ∧ E where L is a conjunction of linear equations and E a conjunction of product-zero-equations. Elementary Flux Modes. The support of a variable assignment α :

V → R with V ⊆ V is supp(α) = {x ∈ dom(α) | α(x) = 0}
. Given a linear system φ, the EFMs of φ are the minimal support solutions of φ over R + . The R-EFMs of φ are the minimal support solutions of φ over R. Note that the interpretation of R is natural for the steady-state equations of metabolic networks with reversible reactions, while the reactions of our networks are always irreversible.

Theorem 20 [START_REF] Allart | Computing sign abstractions of linear systems[END_REF]. Let φ be a mixed linear system. We can compute in at most exponential time some conjunctive formula φ with existential quantifiers such that h B • sol R+ (φ) = sol B (φ ).

Proof sketch. There are quite some insights behind this theorem that we can only sketch here. First, any linear equation L system can be rewritten in the form Ay . =0 where A is an integer matrix and y a sequence of pairwise distinct variables such that V (y) = V (L). Let P be a positive integer matrix whose columns contain all the EFMs of A. These can be computed in at most exponential time [START_REF] Zanghellini | Elementary flux modes in a nutshell: Properties, calculation and applications[END_REF]. Then sol R+ (Ay . =0) is equal to sol R+ (∃x. P x . =y) given that any solution of Ay . =0 over R + can be obtained from some linear combination of the EFMs of A. Second, a formula φ ∈ F Σ is called h B -exact if sol B (φ) = h B • sol R+ (φ). Unfortunately, not every linear systems is h B -exact. However, the formula φ equal to ∃x.P x . =y can be shown to be h B -exact, roughly since matrix P contains only positive integers. Third, it was noticed that any conjunction of product-zero equations is h B -exact as well. Fourth, for any system of product-zero equations E and any sequence of variables z, the formula φ equal to ∃z.φ ∧ E can be shown to be h B -exact. Finally, any mixed linear systems φ can be brought into the form of φ by computing the EFMs of the matrix A of the linear subsystem of φ in exponential time.

Exact Algorithm.. In order to compute the h B -abstraction of a mixed linear system φ, we first compute φ along the lines of the sketch of the proof ideas of Theorem 20. Second, given that φ is a conjunctive formula, we compute sol B (φ ) by finite domain constraint programming.

Characterizing ∆ 3 -Abstractions. We present a characterization of ∆ 3 -abstractions of linear equation systems and show that it provides an exact algorithm solving the objective in the case of ∆ 3 . We first decompose h ∆3 into h B and the binary relation, that is defined by the following first-order formula in the logic with pairs F We next define applications of proj G in FO-logic. For any sequence of variables y and FOformula φ(y) ∈ F 2 Σ with V(φ(y)) ⊆ {y} we define a formula proj G (φ(y)) ∈ F 2 Σ describing the application of proj G to the solutions of φ(y) by ∃z. φ(z) ∧ proj G (z, y) where φ(z) is obtained from φ(y) by replacing the variables in y by arbitrarily but fixed fresh variables z. Theorem 23. For any linear formula L ∈ F Σ we can compute in at most exponential time a positive conjunctive formula with existential quantifiers φ ∈ F Σ such that:

h ∆3 • diff(sol R+ (L)) = ν -1 (sol B (φ))
Proof Let L(y) be a linear system with V(L(y)) = V(y) where y is a sequence of distinct variables. The time for computing φ is dominated by the time for computing the elementary Fig. 8. Predictions for the networks analysed in this paper, obtained respectively by pure abstract interpretation, the heuristic based on minimal support consequences and the exact algorithm.

However, the heuristics with minimal support consequences is not always exact: we found a counter example given on the right for which it slightly overapproximates the exact solution set.

On the other hand, the heurstic algorithm with EFMconsequences is remarkably faster than the exact algorithm -in the benchmark on leucine overproduction, we have 5 minutes versus 5 hours -while still being equally precise in most cases.

First-order expressions

  and formulas: e ∈ EΣ ::= x | c | e e where ∈ F (2) , c ∈ C φ ∈ F Σ ::= e . =e | ∃x.φ | φ ∧ φ | ¬φ where x ∈ V Interpretation of expressions as sets of elements e S,α ⊆ dom(S), where S is a Σstructures and α : V → dom(S) where V contains all free variables. c S,α = c S x S,α = {α(x)} e e S,α = ∪{s S s | s ∈ e S,α , s ∈ e S,α } Interpretation of formulas as truth values φ S,α ∈ B:

Fig. 7 .

 7 Fig.7. Inputs of our algorithm on the benchmark example from of leucine overproduction.

2 +

 2 interested in interpreting this formula over R 2 + .Lemma 21. The relation projG R is a function satisfying h ∆3 = h 2 B • proj G R 2 + .

Lemma 22. proj G R 2 +

 2 • 2-sol R+ (φ(y)) = 2-sol R+ (proj G (φ(y)))

modes, which can be done in at most exponential time.

where φ is an equivalent conjunctive formula for the mixed linear sytem ν(∃z. L 2 (z) ∧ proj G (z, y))

Note that sol B (φ) can be computed by finite domain constraint programming. This yields an exact algorithm for computing the ∆ 3 -abstraction of a system of linear equations, which is a special case of our general objective without kinetic constraints.

For adding a treatment of kinetic constraints, we consider the union B ∪ ∆ 3 as a relational structure providing the values and functions of both structures B and ∆ 3 . The signature of this mixed structure consists of the function symbols in {+ B , * B , + ∆3 , * ∆3 } and the constants in the set B ∪ ∆ 3 , all of which are interpreted by themselves in the mixed structure B ∪ ∆ 3 . The set of first-order formulas over the mixed signature is denoted by F B∪∆3 . For any α : V → dom(S), we can define its restriction α |V : V → dom(S) such that for all y ∈ V ⊆ V , α |V (y) = α(y).

Proposition 24. For any formulas φ ∈ F Σ and φ ∈ F Σ∪∆3 and sets of variables V ⊆ V = V(φ) ∪ V(φ ) we can compute in linear time a formula φ M ∈ F B∪∆3 over the mixed signature such that:

The set sol B∪∆3 (φ M ) can be computed by a finite domain constraint programming, since B∪∆ 3 is a finite structure. By combining Theorem 23 and Proposition 29 we obtain an algorithm for solving the general problem of Section 8 in the cases of ∆ 3 .

Characterizing ∆ 6 -Abstractions. The case of ∆ 6 is considerably more envolved that the case of ∆ 3 , even though follwoing the same general approach. For this, we consider the abstraction h ∆6 as an element of the algebra of total functions on R 2 + , that we denote as R 2 + → R 2 + . The following lemma shows that h ∆6 is the sum of h ∆3 and h 2 B in this Σ-algebra.

). Furthermore, we define for any two functions g :

We can now define the ternary relation idproj G R 2

+ in the first-order logic with pairs by idproj G :

Σ such that for all x, y 1 , y 2 ∈ V:

We next define applications of idproj G in FO-logic. For any sequence of variables y and FO-formula φ(y) ∈ F 2 Σ with V(φ(y)) ⊆ {y} we define a formula idproj G (φ(y)) ∈ F 2 Σ for describing the application of idproj G to the solution set of φ(y). We let idproj G (φ(y) be ∃z. φ(z) ∧ idproj G (z, y 1 , y 2 ) where φ(z) is obtained from φ(y) by replacing the variables in y by arbitrarily but fixed fresh variables z and by fixing two sequences of fresh variables y 1 , y 2 ∈ V m such that y = (y 1 , y 2 ).

Theorem 28. For any linear formula L(y) with free distinct variable y we can compute in at most exponential time a positive conjunctive formula with existential quantifiers φ ∈ F Σ and sequences of variables y 1 , y 2 such that:

Proof Let L(y) be a linear formula L with free distinct variable y ∈ V m .

We can now finish the proof by computing the h 2 B abstraction of the above solution set similarly to the case of ∆ 3 .

where φ is conjunctive formula equivalent to the mixed linear system ν(∃z. L 2 (z) ∧ idproj G (z, y))

For adding a treatment of kinetic constraints, we consider the union B ∪ ∆ 6 as a relational structure providing the values and functions of both structures B and ∆ 6 in analogy to the case of ∆ 3 . We denote the of first-order formulas over the signature of this mixed structure by F B∪∆6 .

Proposition 29. For any formula φ(y) ∈ F Σ with distinct free variables y, formula φ ∈ F Σ∪∆6 , and variable sets

The set sol B∪∆6 (φ M ) can be computed by a finite domain constraint programming, since B∪∆ 6 is a finite structure. By combining Theorem 28 and Proposition 29 we obtain an algorithm for solving the general problem of Section 8 in the case of ∆ 6 .

Heuristic Algorithm Based on Minimal Support Consequences

The intuition behind the heuristic with minimal support consequences relies on two facts: first, adding consequences to a given linear system L before applying abstract interpretation can improve the precision of the abstraction, as already discussed (2); second, the smaller the number of variables in an equation, the more constraining generally is its abstract interpretation. The heuristics is therefore very simple: before abstracting from R + to ∆, the linear system L containing the steady state equations of the system is replaced by a linear system L min containing all the minimal support R + -consequences of the equations in L. The linear system L min can be computed by applying any existing algorithm for the calculation of R-EFMs to the orthogonal complement of L ⊥ as follows:

1. From L compute a linear system L ⊥ whose solution space -seen as a subspace of the vector space R V(L) -is the orthogonal complement of sol R+ (L). This can be done for example by using a variant of Gauß' triangularization method. 2. From L ⊥ compute the R-EFMs l min 1 , . . . , l min k with any known R-EFMs algorithm. 3. Build L min by using l min 1 , . . . , l min k as the coefficients of the equations of L min .

Experimental Results

We experimentally compare three algorithms for overapproximating the objective {β |V | β ∈ sol R+ V (φ) ∆6 ∩sol ∆ V (φ )} given a linear system of equations φ, kinetic constraints φ , and observable variables V . The first algorithm directly applies pure abstract interpretation to φ to compute sol ∆6 (∃V .(φ ∧ φ )) by finite domain programming where V = V(φ) \ V , overapproximating the objective by Theorem 18. The second algorithm enriches the linear system φ with its minimal support consequences as discussed in Section 10 before applying abstract interpretation. The third algorithm is the exact algorithm that can be derived from Theorem 28.

The experimental results are summarized in Fig. 8. The first instance verifies our expectations on the toy metabolic network with a simple cycle in Fig. 1, without kinetic constraints and V = {X, Y } as observable variables. The exact algorithm shows that there are 6 abstract solutions, one for each value of ∆ 6 . The heuristic with minimal support consequences finds exactly these same 6 abstract solutions, while by applying pure abstract interpretation we find 19 abstract solutions (out of 36 possible assignments), thus a large overapproximation.

The second real scale instance treats leucine overproduction on the network from Fig. 2, see Fig. 7 for a discussion of the precise inputs. The heuristic and the exact algorithm produce the same result with 226 abstract solutions, while by pure abstract interpretation 292 abstract solutions are found. Thereby, the two new algorithms both remove the same 2 wrong gene knockout predictions with respect to baseline algorithm by pure abstract interpretation.