William Briand

Marc Rébillat

Mikhail Guskov

Nazih Mechbal

Damage imaging post processing for delamination size assessment of CFRP aeronautic structures

Keywords: Structural Health Monitoring, Quantification, Composite, Guided waves, Piezoelectric

Thanks to their high strength to mass ratio, composite materials are now widespread in the aerospace industry. Nevertheless, this type of material is subject to internal damages like delamination. In order to detect and localize these damages, robust and precise Structural Health Monitoring algorithms exist for this purpose and have been validated experimentally. However, in order to avoid structures catastrophic failures and to estimate their residual life, there is still a huge need of reliable damage size assessment methods. In this paper, a damage quantification method is proposed. This strategy is based on the extraction of a damage size sensitive feature computed from damage imaging results. Here damage imaging stands for methods that use ultrasonic Lamb waves-based map of damage localization likelihood index. This feature is extracted from each labelled example of a training set in order to infer a mathematical model used to predict the area of a delamination of unknown damages. The proposed method is successfully validated on experimental data carried out on CFRP plate samples equipped with a piezoelectric transducers network.

Introduction

In order to reduce their maintenance costs, airlines are increasingly interested in predictive maintenance systems. This type of maintenance consists of immobilizing an aircraft once a condition indicator exceeds a threshold, rather than wait for a pre-determined period of use. This indicator reflects the actual state of degradation of the monitored structure. To build this indicator, it is necessary to know the size of the detected damages to be able to estimate the remaining life of the structure. This is why there is a strong need for reliable and robust quantification algorithms in SHM field. In this paper, a new quantification method based on post-processing of images resulting of localization algorithms is proposed. The localization methods of interest here are the ones based on Lamb waves propagation by means of a piezoelectric elements network. Some attempts have already been carried out to post-process Lamb wave based imaging results for damage size quantification purposes. In [START_REF] Kulakovskyi | Development of a SHM System by Elastic Guided Waves Applied to Aeronautic Structures[END_REF], the authors proposed a data-driven approach using a Convolutional Neural Networks trained on a dataset composed of images generated by SFE method and a localization algorithm called Excitelet [START_REF] Quaegebeur | Dispersion-based imaging for structural health monitoring using sparse and compact arrays[END_REF]. In this dataset, an aluminium plate contains a hole with various sizes and positions. Once the model performs well on this training set, it is applied on numerical and experimental unknown datasets. Another quantification approach aims to assess the size of a crack and a hole [START_REF] Migot | Numerical and experimental investigation of damage severity estimation using Lamb wave-based imaging methods[END_REF]. It consists in localizing the tips of a crack and thus assessing its length, or localizing the edge of the hole and estimating its diameter using two different imaging techniques. The method is applied on an aluminium plate. In [START_REF] Sorrentino | Improved elliptical triangulation method for damage detection in composite material structures[END_REF], the authors developed an algorithm based on Time of Arrival (ToA) localization to assess the size of an impact [START_REF] Sorrentino | Improved elliptical triangulation method for damage detection in composite material structures[END_REF]. Each tip of the damage is localized and the damage size is computed as the area of the polygon formed by these tips. An application have been made on composite plate with data from simulation and from experiments. However, in these articles only one localization method is applied. Thus the applications are done on simple aluminium plate with artificial damage and rarely on CFRP material with delamination. Moreover, there is no method in the literature that can post-process images from different localization techniques and compare the results with each other. The approach proposed in this paper consists in monitoring the growth of a damage based on the computation of a new damage index that varies in a monotonic manner with the size of the damage. This feature is extracted from images coming from damage localization methods: it corresponds to the area of the region surrounding the estimated position of damage where the damage index is higher than a threshold. This damage index is computed for each element of a training set composed of previous damaged states of the structure. Then a data-driven model is built with a polynomial regression. Finally, the size of unknown damages is assessed. Results of this method are compared for four different damage localization techniques: Time of Arrival (ToA) and Difference Time of Arrival (TDoA) [START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF], Delay And Sum (DAS) [START_REF] Michaels | Detection, localization and characterization of damage in plates with anin situarray of spatially distributed ultrasonic sensors[END_REF] and Reconstruction Algorithm for the Probabilistic Inspection of Damage (RAPID) [10]. As more and more aeronautic structures are made of composite materials, a reliable quantification method must be tested on this type of material. The proposed method is successfully validated on experimental data coming from CFRP plate samples equipped with a piezoelectric transducers network and undergoing fatigue solicitations. A method to automatically select its own parameters is proposed. Prediction results with a mean error of 6% are obtained on experimental results with the proposed approach.

In the following, the feature used for for the post-processing method is explained in details. Then proposed approach is explained in details. Finally, the post-processing method is successfully validated on experimental data of fatigue tests from NASA prognostics data repository with delamination close to real life damage.

Theory

HDLI feature

The new damage index proposed in this paper is based on post-processing of localization results. Typically, a localization method takes raw signal from piezoelectric sensors as input, and process them to obtain a map of the structure where each point is associated with a value of a particular Damage Localization Index (DLI). The higher this index is, the higher is the likelihood of finding a damage at this position. In this article, we will compare the performance between different common localization algorithms: Time of Arrival (ToA) and Difference Time of Arrival (TDoA) [START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF], Delay And Sum (DAS) [START_REF] Michaels | Detection, localization and characterization of damage in plates with anin situarray of spatially distributed ultrasonic sensors[END_REF] and Reconstruction Algorithm for the Probabilistic Inspection of Damage (RAPID) [10].

The idea is to compute a feature based on the previously mentioned maps that varies in a monotonic manner with the damage size. Thus the damage index is defined as this area surrounding the estimated damage position where the DLI is higher than a threshold. The assumption made here is that the area of this region depends on the actual size of the damage. To perform the computation of the area considered, a threshold is applied to the map giving a binary image where only the region with a DLI higher than the threshold can be identified. Then an image segmentation algorithm selects the region surrounding the estimate damage position and computes the area A of the region. Finally, the feature is defined as the ratio between the area calculated and the total area of the structure A tot . This damage index is called High Damage Localization Index (HDLI):

HDLI = A A tot (1)

Damage quantification model

The proposed approach aims to estimate the size of an unknown damage. To perform this task, a quantification model is built based on the HDLI feature previously described. This model is data-driven thus the workflow is divided into a learning step and a prediction step. The learning set is composed of N damage cases labeled with the associated damage size {(S 1 , s 1 ), . . . , (S N , s N )}. S i is the set of signals from the transducers for the i th damage case and s i is the corresponding size of the damage.

Then each example of this dataset is processed with a given localization method, resulting a new dataset {(I 1 , s 1 ), . . . , (I N , s N )} where I i is the image returned by the localization algorithm. The HDLI feature h i is computed for each image, resulting in the learning dataset X = {(h 0 , s 0 ), . . . , (h N , s N )} that will be used to train the regression model. HDLI values close to 0 are discarded to improve the sensitivity of the method for large damages since the purpose is to extrapolate a model toward higher damage sizes. In the following, H ∈ R N and S ∈ R N will denote respectively the vector of HDLI features and the vector of damage sizes of all N training samples. A polynomial regression of order d is then performed on X.

Ŝ = d j=0 β j H j (2)
In Eq. 2, d is the degree of the polynomial, Ŝ ∈ R N is the vector of estimated damage size and H j ∈ R N is the vector of j th power of the components of H. To avoid over fitting of the regression problem, a variation of classic linear regression called ridge regression is used here [START_REF] Friedman | The Elements of Statistical Learning[END_REF]. It consists in adding a penalty term λ on the parameters β j in the ordinary least square regression problem in order to control for their amplitude. The optimal λ value is chosen to minimize the error over the learning set using a gradient descent algorithm. To avoid overtraining of the model on the learning data, a crossvalidation technique called Leave-One-Out is used. One must chose a value for the degree of the polynomial. Any degree high enough will yield to the same regression model because coefficients relative to high degree terms will be close to zero thanks to the regularization. Throughout this paper, d is set to 8.

The data-driven model built earlier is then used to predict the size of an unknown damage. The first step is to process the signal with the damage localization method to get a DLI image of the structure. Then the HDLI h test is computed from this image. The corresponding size s test is finally estimated using the size quantification model previously inferred.

Methodology

The overview of the method's steps are:

Step 1 Get signals from each transducer on the plate.

Step 2 Process the signals corresponding to the i th case with one of the damage localization algorithm described earlier. The result is a DLI image of the structure.

Step 3 Compute the HDLI associated with the image i.

Step 4 Repeat steps 2 and 3 for each i ∈ [[1; N ]].

Step 5 Perform the polynomial regression model using the training set of HDLI previously computed and the known damage size of the corresponding damage case.

Step 6 Compute the HDLI value of an unknown damage case.

Step 7 Use the regression model to estimate the size of the damage.

In order to assess the performance of the inferred model, the following metric will be used in the rest of this paper. It is defined as the a relative error computed on learning or test set

ε = 1 n L i=1 |ŝ i -s i | s i ( 3
)
where L is the size of the dataset (learning or test set).

Parameters tuning

Table 1: Localization parameters selected for the rest of the paper.

Loc. method Method parameter Symbol Value

ToA & TDoA Decay rate of an exponential windowed function applied to reduce secondary reflections. The model described earlier depends on a set of parameters. In order to properly quantify the size of a damage, the data-driven model need to associate a HDLI value to a unique damage size. Thus the tuning of the model parameters need to result in a monotonic manner of the model. As mentioned earlier, 4 different localization methods are tested in this paper. Each of this algorithm depends on a tuning parameter. The authors have found that these parameters have little impact on the quantification model, thus the value are selected empirically in Table 1.

The post-processing method presented in this paper aso depends on a threshold level called T in the HDLI computation step. The optimal threshold value T opt is selected as the one minimizing the error of the training set:

T opt = argmin T ∈R 1 2N N i=1 s i -ŝi 2 (4)
Since this cost function is not smooth, the derivative free optimization method called simplex [START_REF] Lagarias | Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions[END_REF] is used.

Application to experimental data

The post-processing method described in the previous section is applied to experimental data. These data are provided by NASA and Standford University [START_REF] Saxena | CFRP Composites Data Set[END_REF].

Data description

The experimental setup consists in a CFRP sample with a dog-bone shape notched at the narrow part to induce a delamination in this region. 12 PZTs are bonded to the sample (6 on each tap). The PZTs placed at the top act only as actuators whereas the ones at the bottom acts only as sensors. The sample is placed in a tensile machine and tested in fatigue at a frequency of 5Hz and a stress ratio of 0.14. The experiment is regularly interrupted to perform Lamb waves propagation measurements along with X-ray control. From theses X-ray images, the delamination can be seen growing progressively with the number cycles. The area of the delamination is extracted from these images with the software Digimizer 1 . This leads to a dataset composed of a collection of raw signals labeled with the corresponding size of the damage. In the following, the first 70% of the dataset is used as learning set and the last 30% is used as a test set to measure the performance of the model on unknown damages.

The samples under study are the L1 coupons from the NASA database, corresponding to the stacking sequence [0 2 /90 4 ] s . The excitation signal of the PZT is a 5-cycles tone burst with a central frequency of 250kHz and an amplitude of 50V. This frequency is selected as it's the one where the Lamb waves modes S0 and A0 are the most distinct as mentioned by [START_REF] Larrosa | In situ damage classification for composite laminates using Gaussian discriminant analysis[END_REF].

Results and analysis

The image post-processing quantification method described in this paper is applied with the four different localization algorithms mentioned earlier. For each sample and each localization algorithm, the performance of the quantification process is evaluated by two means. First, a graph with the true size of the delaminated area on X axis and the size estimated with the quantification process on Y axis is shown Fig. 1. For the sake of space, L1 S 18 results are not shown under graph form since the scale is different from that of the other coupons. The X = Y line corresponds to a prediction without error. Thus the closer the points are to this line, the more accurate is the prediction. The other way to measure the performance of the algorithm is to compute the error 3 on both learning and test set. ToA and TDoA exhibit poor results on both learning and test set. DAS method shows good performance on learning set but mediocre results for prediction on the test set. Finally, the RAPID method performs well on both learning and test sets. This tendency are confirmed both the error obtained table 2. It seems clear that the selected localization method highly influence the results of the quantification process. The HDLI feature computed with ToA and TDoA show no correlation with the true damage size, leading to poor performance on both datasets. An explanation could be that this localization methods only deal with the time of flight which is not very influenced by the size of the damage. On top of that, the high DLI region is made up of ellipses or hyperbolas which area does not clearly vary with the damage size. This phenomenon could influence the performance of the regression as well. Besides, DAS and RAPID both exhibit great results on the training set. This performance could be explained by the fact that DAS and RAPID take amplitude of the signals into account in addition to time of flight, which increase the dependence of the HDLI feature to the actual size of the damage. Moreover, high DLI area in the DAS and RAPID methods is circular and HDLI varies in a monotonic manner over the size of the damage which lead to a more accurate model. DAS method shows poor results on testing set. Indeed, above a certain level, the HDLI remains steady with the damage size. It means that above this level, each damage size leads to the same HDLI. On the opposite, RAPID method leads to a monotonic function between HDLI and damage size for the full range of learning set and test set.t

Conclusion

In this paper a new quantification method for damage size assessment of PZTs instrumented structure was presented. This new algorithm consists in post-processing the image results from a given localization method to compute a feature called HDLI. This damage index is correlated with the actual size of the damage. A quantification model is then inferred on a learning dataset. The model thus obtained is used to predict the size of unknown damage. This method was tested with experimental data coming from fatigue tests of CFRP samples. The proposed approach was tested on four common localization methods and show promising results with a mean error of 6% obtained. In the current approach, it is necessary to train a single model for each coupon under testing. Future work will focus on transfer learning, i.e. learning the model on the full dataset available for a given coupon and performing prediction with the quantification model on an another sample with similar material and geometrical properties.

τ 5 .0 10 - 6 DAS Number of samples over which time integration is performed. K 1 RAPID

 561 Parameter set to adjust the spread of the ellipses around each path. β 1.05

Fig. 1 :

 1 Fig. 1: Damage quantification performance using different damage localization methods on L1 coupons data from NASA. ToA (a), TDoA (b), DAS (c), RAPID (d). The diagonal Y = X line correspond to a prediction without error.

Table 2 :

 2 Performance over the dataset measured with training and test error.

	Dataset ToA	TDoA	DAS	RAPID
	ε train εtest ε train εtest ε train εtest ε train εtest
	L1 S11 200 32 140 36 36 31 12 7.8
	L1 S12 220 50 300 70 62 38 13 6.2
	L1 S18 390 49 410 36 31 20 19 7.8
	L1 S19 750 85 400 78 29 47 38 27

https://www.digimizer.com/

Funding

This work has received funding from the European Union's Horizon 2020 research and innovation program under the REMAP project (grant agreement number 769288). https://h2020-remap.eu/


Joachim Niehren

Cristian Versari

Mathias John

François Coutte

Philippe Jacques

Predicting Changes of Reaction Networks with Partial Kinetic Information

Keywords: Reaction networks, model based prediction, abstract interpretation, constraint solving, metabolic engineering, genetic engineering

come   L'archive ouverte pluridisciplinaire

Introduction

Mathematical methods for analysing reaction networks [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors-I. the deficiency zero and deficiency one theorems[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] often require full kinetic information for all reactions, while in systems biology practice only partial information is usually available. Therefore, we study the problem of how to model reaction networks with partial kinetic information, and how to reason qualitatively about such models with methods from computer science. In particular, we wish to predict changes of reaction networks with partial kinetic information that may or must lead to a target change of their steady state.

When full kinetic information is not available, the existing model-based reasoning methods tend to ignore the kinetic information all over. Most typically, this holds for 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] when applied to metabolic networks [START_REF] Otero | Industrial systems biology[END_REF][START_REF] Sohn | In silico genome-scale metabolic analysis of Pseudomonas putida kt2440 for polyhydroxyalkanoate synthesis, degradation of aromatics and anaerobic survival[END_REF]. The missing information is then compensated heuristically by the adoption of ad hoc optimization criteria. Alternatively, pathway analysis approaches [START_REF] Papin | Comparison of network-based pathway analysis methods[END_REF] rely on the structure of reactions networks, but the combinatorial nature of the problem makes difficult their application to densely interconnected networks. Both above methods have extensions that deal with partial kinetic information about inhibitors. This is done by adding boolean constraints that state the conditions on which an inhibitor blocks a reaction [START_REF] Jungreuthmayer | Designing optimal cell factories: integer programming couples elementary mode analysis with regulation[END_REF]. However, blocking inhibitors are not really appropriate in deterministic semantics, where the average over blocked and unblocked situations is to be considered. Therefore, it remains open how to deal with nonblocking inhibitors, which only slow down a reaction in average.

In the first part of this article, we propose a new modeling language for reaction networks with partial kinetic information, a short version of which was presented at the CMSB'2015 conference [START_REF] Niehren | Qualitative Reasoning about Reaction Networks with Partial Kinetic Information[END_REF]. Our language is parameterized by a similarity relation on kinetic functions, so that the rate laws of chemical reactions need only to be specified up to similarity. For instance, two kinetic functions could be considered similar if they have the same monotonic behavior. Then, the kinetic function mapping x to 42x or to 5x/(7 + x) would be similar, since whenever x increases then the values of both terms increase, and whenever x decreases then they both decrease.

Reaction networks are applied to a chemical solution, which is typically placed within a context. This may be another chemical solution that is subject to an adjacent reaction network or to an experimental environment. The situation is illustrated in Figure 1. Since we do not want to model the context, we equip reaction networks with an interface to possible contexts only. The interface specifies which species can inflow from the context and which species may We also assume that some of the reactions of the network may be candidates for knockout. The choice of whether a reaction is knocked out or not remains external to the network. For the prediction task we are interested in network changes possibly combining multiple reaction knockouts and influx changes.

The models of reaction networks in our language have a graphical syntax that is reminiscent of Petri nets, and also an equivalent Xml syntax. They are given a steady state semantics in terms of arithmetic equations, which define the relation between influxes and outfluxes of the network in steady states, depending on which knockout candidates were knocked out. The steady state semantics subsumes the usual flux balance equations, but enhanced with equations on the rates of reactions based on the partial kinetic information. This information is expressed by variables for kinetic functions that are subject to similarity constraints. In this way, the inhibitors of a reaction slow its rate down rather than shutting it down completely, while the activators of a reaction speed its rate up.

Our language can be used to model metabolic networks with complex regulation such as for B. subtilis in the Subtiwiki [START_REF] Mäder | Subtiwiki 2014a comprehensive community resource for the model organism bacillus subtilis[END_REF]. An example is the regulation network of the PIlv-Leu promoter of B. subtilis, which regulates the metabolism of the branched-chain amino acids Valine, Leucine, and Isoleucine. Previous models of these metabolic networks in the Subtiwiki were not given any formal semantics, so that they could not be used directly for qualitative prediction algorithms. A detailed model of this network in the language presented here was published in [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] recently.

In the second part of this article, we develop a new prediction algorithm that can be applied to any reaction network modeled in our language. This algorithm does qualitative reasoning [14] based on abstract interpretation [START_REF] Cousot | Systematic design of program analysis frameworks[END_REF], by which the partial kinetic knowledge is discretized while the unknowns are abstracted away. A first version of the algorithm was presented at the VMCAI'2013 conference [START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF], but restricted to networks of reactions with mass action kinetics. The second version, presented at the CMSB'2015 conference [START_REF] Niehren | Qualitative Reasoning about Reaction Networks with Partial Kinetic Information[END_REF] and extended here, removes this limitation.

Given a reaction network, abstract interpretation can be applied to infer a difference constraint that relates network changes to changes in the steady states. This can be used to predict which network changes may or must lead to an expected change of the outfluxes. As already stated above, the network changes that we are interested in are reaction knockouts, influx changes or multiple combinations thereof. This becomes possible since the models of reaction networks in our language do describe how network changes affect the steady state semantics.

Difference constraints are finite domain constraints from a novel class. The second step of our prediction algorithm consists in computing the set of all solutions of a difference constraint. Since difference constraints are finite domain constraints, their solution set is always finite. We built two constraint solvers for difference constraints based on finite domain constraint programming. The solver reported earlier in the conference paper preceding this article was developed from scratch in the programming language Scala [START_REF] Odersky | The evolution of scala: Ple'14 keynote[END_REF]. Since then we developed a new solver by reduction to the MiniZinc solver [START_REF] Nethercote | Minizinc: Towards a standard CP modelling language[END_REF] for finite domain constraints. While our Scala solver could enumerate only the n-best solutions where n ≤ 3000 while consuming considerable time (more than 10 minutes), the MiniZinc solver can indeed return the complete set of all solutions sets in all our applications, while enumerating more than 5000 solutions in less then a second.

Solutions of difference constraints can be interpreted as predictions of network changes that may or must lead to an overproduction target. When only seeing the n-best solutions, one can find solutions with few network changes that are compatible with the overproduction target. But since we are now having access to the complete solutions sets of difference constraints obtained from reaction networks in practice, we can distinguish the solutions that are merely compatible with the overproduction target from those that safely entail it. As we will argue, safe solutions correspond to predictions that must satisfy the overproduction target, while compatible solutions yield predictions that may satisfy the overproduction target or not. It also turns out that multiple network changes may be necessary to obtain safe solutions, while single knockouts may not be enough.

We have implemented the prediction algorithm and integrated it into a prediction tool. We illustrate this tool and our prediction algorithm at the example of two simplified models of leucine production of B. subtilis, that focus on the regulation of the PIlv-Leu promoter. The target here is leucine overproduction. For the simpler of the two models, the prediction based on the graphical model can be done manually by humans. This illustrates that our algorithm formalizes a natural kind of qualitative reasoning.

In a follow up work [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF], our algorithms were applied to single knockout prediction for leucine overproduction in a larger and more realistic model. The predictions obtained there were not safe, but still 6 out of 14 predictions could be verified successfully by gene knockouts in the wet lab.

Compared to the previous two conference papers [START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF][START_REF] Niehren | Qualitative Reasoning about Reaction Networks with Partial Kinetic Information[END_REF], we extended the prediction algorithm to multiple changes, rather than single reaction knockouts or single influx changes. We also introduce the notion of safe predictions, for which the complete set of solutions of a difference constraint must be considered. While safe solutions could have been defined in theory before, their practical relevance became apparent only with the new MiniZinc constraint solver for difference constraints, which provides very efficiently the enumeration of the set of all solutions of difference constraints inferred from reaction networks containing up to about a hundred reactions.

Part I

Modeling Language

Reaction networks with complete kinetic information are introduced in Section 2, while a language for describing such networks modulo a similarity relation on kinetic functions is introduced in Section 3. In Section 4, we illustrate the language for the modeling of two simplified models of leucine production and informally explain how steady state equations can be exploited for change prediction. In Section 5 we define appropriate similarity relations on kinetics functions to be used in our modeling language.

Reaction Networks

We define reaction networks with complete kinetic information, and show how to compute their steady state semantics. The notion of reaction networks and how to infer their ODEs is basically standard (see e.g. [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors-I. the deficiency zero and deficiency one theorems[END_REF]), may be with the exception of inflows, outflows, and knockout candidates.

We want to model biological systems that can interact with their context, i.e., their biological or experimental environment as illustrated in Figure 1. For instance, the chemical solution can model all species of a cell, and the context the cultivation medium of the cell for feeding or extracting species. Alternatively, the chemical solution can represent the subset of species affected by a selected pathway of a cell, and the context all species that are produced and consumed by connected pathways.

In order to model the interaction of a system with all its possible contexts, we add inflows and outflows to our notion of reaction networks. An inflow increases the concentration of some species in the system, like a reaction that produces a species from nothing. But in contrast to all other reactions, the kinetic rates of inflows are left completely free according to the assumption that they are controlled by the context, which may vary from case to case. Therefore the rates of the inflows are unknown to the system itself. An outflow decreases the concentration of some species of the system. For simplicity, we assume that the kinetic rates of the outflows always follow the mass action law with some parameters fixed by the system. An outflow is thus like a reaction with mass action kinetics that consumes a species while producing nothing. Alternatively, outflows could also be modeled by distinguishing a subset of reactions with fixed kinetic rates.

Furthermore, we also want to reason about a set of possible biological systems at the same time, which differ only in one or many reaction knockouts. Since knockouts do not make sense for arbitrary chemical reactions, each network will distinguish a subset of reactions that are candidates for knockout. In the ODEs of reaction networks, we will then use a Boolean variable for each knockout candidate for selecting whether a candidate is knocked out or not.

We denote the set of non-negative real numbers by R + and the set of Booleans by B = {0, 1}. Note that we assume that B ⊆ R + . A kinetic function of arity k ≥ 0 is a function of type κ : R k + → R + . Kinetic functions will be used to define the rate laws of chemical reactions. Let S be a finite set of species. A chemical solution with species in S is a function of type S → R + that maps each species to its concentration.

A chemical reaction r is a tuple of the form: s 1 , . . . , s k κ -→ s k+1 , . . . , s l where 0 ≤ k ≤ l, s 1 , . . . , s l ∈ S, and κ : R k + → R + is a kinetic function. Any reaction has a tuple of reactants s 1 , . . . , s k and a tuple of products s k+1 , . . . , s l . In order to account for the stoechiometry of a reaction, we write rct r (s) to denote the number of occurrences of s in the tuple of reactants of r, and prd r (s) for the number of occurrences of s in the tuple of products of r. A modifier of a reaction is a species s with rct r (s) = prd r (s). Whether a modifier behaves as an activator or as an inhibitor depends on the choice of the rate law κ. Definition 1. A reaction network is a tuple N = (S, R, I, O, K) where S is a finite set, R is a finite set of chemical reactions with species in S, K ⊆ R a subset of knockout candidates, I ⊆ S a subset of inflow species, and O ⊆ S ×R + a partial function mapping outflow species to their rate constant.

An inflow species is an element of I and an outflow species an element of the domain dom(O). Each inflow species s ∈ I specifies an inflow that increases the concentration of s in the chemical solution. Each outflow species s ∈ dom(O) describes an outflow into the context that decreases the concentration of s in the chemical solution. While the rates of inflows are controlled by the context, the rates of outflow species depend on the concentration of the outflowing species. We assume that the rates of outflows follow the mass action law with constant O(s), i.e. it is equal to O(s) • z s , where z s is the concentration of s in the chemical solution. We note that any species may have an inflow and an outflow at the same time.

For any fixed context fixing the rates of the inflows, a reaction network defines the evolution of the concentration of all species of a chemical solution over time. We omit the formal definition, since we are exclusively considering steady states. Steady states are limit time points where the concentrations of all species, the rates of all reactions, inflows and outflows have become stable. It should be noticed that steady states of reaction networks are unique if they For all species s ∈ S define arithmetic expressions:

rct(s) = df r∈R rct r (s) • v r prd (s) = df r∈R prd r (s) • v r cons(s) = df x s + rct(s) if s ∈ I rct(s) otherwise prod (s) = df y s + prd (s) if s ∈ dom(O) prd (s)
otherwise Inference rules for arithmetic equations: exist, for any given initial solution and any context. However, since we neither fix the initial solution nor the context, many different steady states may exist for the same reaction network.

(rate outflux ) s ∈ dom(O) y s = O(s) • z s (rate ¬ko ) r ∈ R \ K r is equal to s 1 , . . . , s k κ -→ s k+1 , . . . , s l v r = κ(z s1 , . . . , z s k ) (rate ko ) r ∈ K r is equal to s 1 , . . . , s k κ -→ s k+1 , . . . , s l v r = o r • κ(z s1 , . . . , z s k ) ∧ o r ∈ {0, 1} (flux balance) s ∈ S cons(s) = prod (s)
Our next objective is to define the steady state semantics of a reaction network by a system of arithmetic equations. These equations will be built over the following set of variables, for values in R + :

• for any species s ∈ S there is a variable z s that denotes the concentration or activity of s in a steady state;

• for any inflow species s ∈ I there is a variable x s that denotes the influx of s, i.e., the rate of the inflow of s;

• for any outflow species s ∈ dom(O) there is variable y s for the outflux of s, i.e., for the rate of the outflow of s in a steady state.

• for any reaction r ∈ R there is a variable v r for the rate of r in a steady state;

• for any knockout candidate r ∈ K there is a variable o r whose value will be 1 if r is knocked out and 0 otherwise.

We assume that the set of variables is totally ordered, so that any subset of variables can be identified with a tuple. The steady state equations of a reaction network are defined in Figure 2. For any species s ∈ S, we define the arithmetic expression rct(s) that sums the consumption rates of all reactions of which s is a reactant. The full consumption rate cons(s) adds the influx x s to rct(s) if s is an inflow species, and is equal to rct(s) otherwise. The expression prd (s) sums the production rates of all reactions of which s is a product. The full production rate prod (s) adds the outflux y s to prd (s) if s is an outflow species and is equal to prd (s) otherwise.

There are three inference rules for inferring the steady state equations. Each of them can be seen as an implication, whose conditions are written above the line and whose conclusions are written below. The inference rule (rate outflux ) creates the equation y s = O(s) • z s for any outflow species s ∈ dom(O), i.e., we assume that all outflows satisfy the mass action law. Note that the variables x s for inflow species s ∈ I are completely unconstrained. The rules (rate ¬ko ) and (rate ko ) provide the rate of reaction r by applying its kinetic function to the concentrations of all its reactants. Compared to (rate ¬ko ), the rate in rule (rate ko ) is multiplied by variable o r . The value of this variable is constrained to be Boolean, depending on whether r is knocked out or not. Rule (flux balance) states that consumption and production rates must be balanced for all species in steady states. , which is obtained by projecting the solutions of the steady state equations for N to the following tuple of variables:

• the variables x s for the n inflow species s ∈ I,

• the variables y s for the m outflow species s ∈ dom(O),

• the variables o r for the o knockout candidates r ∈ K.

In the tuple, the inflow variables appear before the outflow variables and these before the knockout variables. The variables of the same type appear in this tuple in the same total order that we assumed on the set of variables.

The variables for concentrations (the z s where s ∈ S) and for rates (the v r where r ∈ R) are implicitly existentially quantified, since they are projected away when defining the exchange relation.

Modeling Language modulo Similarity

We now present a modeling language for reaction networks with partial kinetic information. As a parameter of our language, we assume a similarity relation ∼ on kinetic functions. Rather than specifying rate laws of chemical reactions by kinetic functions, we will describe them only up to similarity: a rate law belongs to ∼κ if it is similar to the kinetic function κ. Enriched chemical reactions will be used to describe the chemical reactions of a reaction network. An example is given in Figure 3. The graph there represents an enriched chemical reaction r with substrate S twice, activator A, an accelerator A , and inhibitor I, and product P once. Please note that the same species may play different roles even in the same reaction, and several times. Activators are like enzymes. Both activators and accelerators speed up a reaction, with the difference that all activators of a reaction must be present for its application, while accelerators do not need to be there.

For graphical representation, we use conventions similarly to Petri nets. Species are represented by rounded nodes s containing the name s of the species, and enriched reactions are graphically represented by boxed nodes r containing the name r of the reaction. Reactions that belong to the set of knockout candidates K are colored in light orange r .

More generally, enriched chemical reactions have different kinds of reactants, that are fixed by a finite set of roles Rol, which is the second parameter of our language. In our example, there will be substrates -that are consumed -and three kinds of modifiers: inhibitors, activators, and accelerators, so we set Rol = {inh, subs, act, acc}. For our graphical syntax, we assign to each role an edge type, for edges pointing from the reactant to the reaction. We will use for subs, for inh, for act, and for acc. The products of a reaction -besides the above modifiers -will be linked by arrows pointing from the reaction to the product. Species may be present n times with the same role. In this case, we annotate the arrows with the multiplicity, as for instance in n . Reactant roles serve to order the arguments of the rate law of an enriched chemical reaction. Such a rate law is given by an enriched kinetic function:

κ : (Rol × R + ) k → R +
We assume that any enriched kinetic function is well-behaved, in that any permutation of arguments with the same role does not change its value. When fixing the order of the arguments, any enriched kinetic function κ can be replaced by a standard kinetic function κ, for instance such that κ(z S , z S , z I , z A , z A ) = κ (subs:z S , subs: z S , inh:z I , act:z A , acc:z A ). An enriched chemical reaction can then be replaced by a chemical reaction, in which the enriched kinetic function is replaced by a variable for a kinetic function. With the same ordering used to get κ from κ , we obtain for the example from Figure 3:

S, S, I, A , A ∼κ --→ P, I, A , A
Here, ∼κ stands for a fresh variable for a standard kinetic function that is similar to κ. A model in our language is a tuple (S, R, I, O, K) where R is a set of enriched reactions, I, O ⊆ S and K ⊆ R. Note that we do not require to specify rate constants for outflows. Graphically, inflow species in I and outflow species in O are indicated respectively by ingoing and outgoing arrows . An example model in graphical syntax is given in Figure 4.

For any model in our language, we can generate a reaction network with variables for kinetic functions that are subject to similarity constraints. Therefore, we can define the steady state equations of any model in the language as before, except that kinetic functions as well as outflows rate constants will be represented by variables. An example is worked out in the next section.

Besides the graphical syntax, our language supports an Xml syntax, which serves for writing the models, so that the graphs can be generated. We implemented tools for doing this in Xslt. These tools can also compute the steady state equations and perform abstract interpretation.

Example: Regulation of Metabolism of B. subtilis

As an example, we model the leucine biosynthesis pathway of B. subtilis in our language. This is one of the complex regulation mechanisms of the metabolism of B. subtilis, for which informal models are given in the Subtiwiki [START_REF] Mäder | Subtiwiki 2014a comprehensive community resource for the model organism bacillus subtilis[END_REF]. The precise similarity relation of the model will be defined in Section 5.

Model Design

Two variations of the model are given in graphical syntax in Figure 4 and Figure 7. The first reaction network describes the base regulation of the ilv-leu promoter (P Ilv-Leu ), and has a single unregulated reaction for Leu production. The second network refines the former, in that a regulation mechanism for Leu production is added, based on the BkL-Bcd operon. Note that we keep the basic model simpler since we do not consider any knockout candidates there, in contrast to the refined model.

For clearer visualization, species nodes have different colors depending on the type of the species: in this paper we will use proteins P , metabolites M , and promoters or binding sites B . The variable z B stands for the activity of the promoter or binding site B, while z P and z M stand for the concentrations of P and M .

We consider an acceleration function with Acc(d) = 1 + d and an inhibition function with Inh(d) = 1/Acc(d). We define the enriched kinetic functions exp such that for all tuples t = (r 1 :

d 1 , . . . , r k : d k ) ∈ (Rol × R + ) k : exp(t) = ri∈{subs,act} d i • rj =acc Acc(d j ) • Inh( r l =inh d l )
Note that the order of arguments with the same role is not important, so that function exp is well-behaved. When a reaction has the exp kinetics, then its inhibitors slow down the reaction but do not block it. Accelerators and activators both speed up the reaction. Furthermore, if one of the activators is missing then the reaction is blocked. One might want to generalize exp with parameters defining the strength of respective accelerations and inhibitions. We do not do so, since these parameters are typically unknown, and since all such generalized expression kinetics will turn out to be similar. Generally, we are only interested in ∼exp, so similar definitions would do the job as well. The enriched mass-action kinetics is the special case ma(t) = exp(t) for all t ∈ ({subs}×R + ) * .

Basic Network

We start with a basic network in Figure 4 which already refines the one from [START_REF] Niehren | Qualitative Reasoning about Reaction Networks with Partial Kinetic Information[END_REF]. 1 Leucine biosynthesis is realized by enzymes which are coded by the genes of the ilv-leu operon. This operon is under the regulation of the promoter P Ilv-Leu . For simplicity, we group the whole reaction network leading to the leucine biosynthesis into reaction r 8 .

The production of leucine depends on the activation of P Ilv-Leu , which is done by reaction r 2 under the regulation of TnrA, CcpA, and CodY. Proteins TnrA and CodY are inflow species added by the context and degraded by reactions r 16 and r 15 respectively. These two proteins play also a regulatory role in leucine degradation by inhibiting reaction r 3 . Protein CcpA is expressed by reaction r 14 and degraded by reaction r 14 .

Transcription at the ilv-leu promoter is well known to be inhibited by CodY through a binding of this latter on the promoter [START_REF] Mäder | Transcriptional organization and posttranscriptional regulation of the Bacillus subtilis branched-chain amino acid biosynthesis genes[END_REF][START_REF] Molle | Additional targets of the Bacillus subtilis global regulator CodY identified by chromatin immunoprecipitation and genome-wide transcript analysis[END_REF][START_REF] Shivers | Activation of the Bacillus subtilis global regulator CodY by direct interaction with branched-chain amino acids[END_REF][START_REF] Tojo | Elaborate transcription regulation of the Bacillus subtilis ilv-leu operon involved in the biosynthesis of branched-chain amino acids through global regulators of CcpA, CodY and TnrA[END_REF][START_REF] Villapakkam | Genetic and Biochemical Analysis of the Interaction of Bacillus subtilis CodY with Branched-Chain Amino Acids[END_REF]. To model this action of CodY on the promoter P Ilv-Leu , we introduce the reaction r 1 which activates the binding site BS CodY of CodY at the promoter, which in turn slows down reaction r 2 and thus reduces the promoter's activity. The binding of CodY to the promoter's binding site BS CodY can be prohibited when CcpA is bound to the promoter. Therefore the presence of CcpA slows down reaction r 1 [START_REF] Brinsmade | Regulation of CodY Activity through Modulation of Intracellular Branched-Chain Amino Acid Pools[END_REF][START_REF] Tojo | Elaborate transcription regulation of the Bacillus subtilis ilv-leu operon involved in the biosynthesis of branched-chain amino acids through global regulators of CcpA, CodY and TnrA[END_REF] but it does not block it on average in a steady state.

The promoter P Ilv-Leu is also down-regulated by Leu in terms of a T-box [START_REF] Brinsmade | Regulation of CodY Activity through Modulation of Intracellular Branched-Chain Amino Acid Pools[END_REF][START_REF] Grandoni | Transcriptional regulation of the ilv-leu operon of Bacillus subtilis[END_REF], which is captured by the negative control of the reaction r 2 by Leu.

Protein TnrA forms a further inhibitor whose impact on the P Ilv-Leu promoter is represented by the binding site BS TnrA through the reaction r 7 , this latter is degraded by reactions r 7 . Protein CcpA is also independently upregulating the ilv-leu operon transcription by binding on it, and thus activating reaction r 2 .The binding of CcpA is captured by BS CcpA through reaction r 9 . This latter is then degraded by r 9 . Moreover, the acceleration of the reaction r 2 by BS CcpA is inhibited when BS TnrA is active. Indeed, the active binding sites BS TnrA and BS CcpA can bind to each other while forming a DNA loop. This phenomenon is captured by the inhibition of the reaction r 9 by BS TnrA .

From the model, the steady state equations in Figure 5 were inferred. These contain variables exp (i) for enriched kinetic functions similar to exp, and variables ma (i) for enriched kinetic functions similar to the mass-action law ma for any i. Note that v r14 = exp (14)

Flux balance equations:

(Leu)

v r8 = y Leu (CcpA) v r14 = v r 14 (CodY) x CodY = v r 15 (TnrA) x TnrA = v r 16 (BS CodY ) v r1 = v r 1 (P Ilv-Leu ) v r2 = v r 2 + v r8 (BS TnrA ) v r7 = v r 7 (BS CcpA ) v r9 = v r 9 Outfluxes:
y Leu = ma (0) (subs: z Leu ) Reaction rates: v r1 = exp (1) (inh:z CcpA , act: z CodY ) v r 1 = ma (1) (subs: (2) (subs:z P Ilv-Leu ) v r7 = exp (7) (act:z TnrA ) v r 7 = ma (7) (subs:z BSTnrA ) v r8 = exp (8) (subs:z P Ilv-Leu ) v r9 = exp (9) (inh: z BSTnrA , act:z CcpA ) v r 9 = ma (9) (subs:z BSCcpA ) v r14 = exp (14) () v r 14 = ma (14) (subs: z CcpA ) v r 15 = ma (15) (subs: z CodY ) v r 16 = ma (16) (subs: z TnrA ) is applied without any arguments. This reflects that r 14 has no reactants or modifiers, so that its rate is constant. Therefore the equations in Figure 5 can be simplified by eliminating the implicitely existentially quantified variables and replacing them by equivalent expressions. In particular, we can simplify the steady state equations in Figure 5 to the equations in Figure 6.

z BS CodY ) v r2 = exp (2) (inh: z BS CodY , acc:z BSCcpA , inh: z Leu , inh:z BSTnrA ) v r 2 = ma

Prediction of Input Changes

In order to illustrate the qualitative reasoning methods that we will develop, we consider the overproduction problem of Leu for the P Ilv-Leu network. The question is: which changes of the influxes may lead to an increase of the Leu outflux? (2) (subs: z P Ilv-Leu ) v r7 = exp (7) (act: z TnrA ) v r7 = ma (7) (subs:z BSTnrA ) y Leu = ma (8) (subs:z P Ilv-Leu ) v r9 = exp (9) (inh: z BSTnrA , act:z CcpA ) v r9 = ma (9) (subs:z BSCcpA ) v r14 = exp (14) () v r14 = ma (14) (subs:z CcpA )

v r2 = v r 2 + y Leu y Leu = ma (0) (subs:z Leu ) v r1 = exp (1) (inh:z CcpA , act: z CodY ) v r1 = ma (1) (z BS CodY ) v r2 = exp (2) (inh:z BS CodY , acc: z BSCcpA , inh: z Leu , inh: z BSTnrA ) v r 2 = ma
x CodY = ma (15) (subs:z CodY )

x TnrA = ma (16) (subs: z TnrA ) Informally, the problem can be solved as follows. Leu is produced only from P Ilv-Leu via reaction r 8 . The rate of this reaction can be increased by increasing the concentration of P Ilv-Leu . Since P Ilv-Leu is solely produced by reaction r 2 , this requires to increase the rate of r 2 . This can be done by either decreasing the concentrations or activities of one of its three inhibitors Leu, BS CodY , or BS TnrA , or by increasing the concentration of its accelerator BS CcpA . Increasing BS CcpA can either be reduced to decreasing the activity of its inhibitor BS TnrA , a choice that we were considering independently already, or by increasing the concentration of its activator CcpA. The latter is not possible, since CcpA is not connected to any inflow, so it cannot be increased by changing the influxes. And the concentration of inhibitor Leu cannot be decreased, since we want to increase its outflux at the same time. Hence, either the activity of BS CodY or BS TnrA must be decreased. This is possible only as follows:

1a. decrease the influx of CodY, or 1b. decrease the influx of TnrA.

Either of these single changes could be enough, but a double change is also compatible with the target of leucine overproduction.

Refined Network

One of the aspects missing in the basic network is that leucine is also degraded for its use in fatty acid biosynthesis. This is regulated by several enzymes captured by the BkL-Bcd operon. This aspect is modeled in the refined network in Figure 7, where we added the species OP BkL-Bcd to the basic network.

We keep the refined network as simple as possible, in order to to illustrate the relevance of multiple change prediction (see [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] for a detailed network). Therefore, we first add some knockout candidates and second add operon OP BkL-Bcd as an inhibitor of the production reaction of Leu r 8 . The activity of operon OP BkL-Bcd is increased by the new reaction r 3 , which is inhibited by TnrA and CodY. Thereby, the presence of both CodY and TnrA now have a negative effect on the production of Leu, while in the basic network there were only positive effects. The second new reaction r 3 degrades the activity of the operon so that it can become steady.

Leu

Prediction of Multiple Knockouts and Input Changes

It turns out that the predictions change with the refined network in Figure 7 by the addition of OP BkL-Bcd .

It is still possible that a decrease of the influxes of CodY or TnrA may lead to an increase of the outflux of Leu, but this must no more be the case. Furthermore, an increase of CodY or TnrA may now also lead to an increase of the outflux of Leu, since it may lead to a decrease of the concentration of OP BkL-Bcd which is an inhibitor of r 8 . This may seem contradictory at a first glance, but in fact the changes of the influxes of CodY or TnrA now have positive and negative effects on the production of Leu. Therefore, none of these changes is safe, in that none of them must necessarily lead to an overproduction of Leu.

When looking into multiple changes, however, which may combine knockouts and influx changes, the contradictory situation for single changes can be resolved. Indeed, it turns out that Leu production will be safely increased by either of the following four double changes: Another way to safely increase the Leu production is to perform a triple change: increase increase but not from zero 3. knockout reactions r 1 , r 3 , and r 7 at the same time.

< = {(x, y) ∈ R 2 + | x < y} decrease > = {(x, y) ∈ R 2 + | x > y} no change . = = {(x, x) | x ∈ R + }
↑ = {(x, y) ∈ R 2 + | 0 < x < y} increase from zero ⇑ = {(0, y) ∈ R 2 + | 0 < y} decrease but not to zero ↓ = {(x, y) ∈ R 2 + | x > y > 0} decrease to zero ⇓ = {(x, 0) ∈ R 2 + | x > 0} no change but not at zero ∼ = {(x, x) | 0 = x ∈ R + } no change at zero ≈ = {(0, 0)}.
The above results will be inferred by our prediction methods for multiple changes in Section 8. There we will also formalize what it means for a change to be safe.

Similarity by Difference Abstraction

We now define similarity relations on kinetic functions by abstracting from changes between real numbers or relations between real numbers. Influx decreases or reaction knockouts are a special case.

We now want to abstract the space of all changes in R + × R + into a finite set of difference relations. For this, we partition the set R + × R + into a finite collection of subsets ∆ ⊆ 2 R+×R+ , so that we can abstract any change in R + ×R + into a unique difference relation of ∆.

In the examples, we will use two different partitions. For predicting inflow changes, we will use the ternary partition ∆ 3 = {<, >, . =}, defined in Figure 8. An influx increase or a speed-up of a reaction can be represented by <, and an influx decrease or a slow-down of a reaction is abstracted to >. The third difference relation . = represents a no change. For predicting reaction knockouts, it is often better to use the senary partition ∆ 6 = {↑, ⇑, ↓, ⇓, ∼, ≈} as defined in Figure 9. Its difference relations can distinguish a reaction slow-down ↓ from a reaction knockout ⇓. Note that ∆ 6 refines the partition ∆ 3 as follows:

< = ↑ ⇑ , > = ↓ ⇓ , . = = ∼ ≈ .
: (R + ) p-1 → R + are similar, written κ 1 ∼ ∆ κ 2 , iff κ ∆ 1 = κ ∆ 2 .
Example 1. Let ma k (subs: d, subs: d ) = k • d • d be the binary mass action law with constant k ∈ R + . For any ∆, the differences abstractions ma k ∆ are equal for all choices of parameter k > 0 since ma ∆3 k = • ∆3 (but not for k = 0). This ternary relation is given by the binary set value function in the table on the right. For instance, the set < • ∆3 > contains <, since an increase (say from 1 to 10) times a smaller decrease (say from 2 to 1) is an increase (from 2 to 20). It also contains > since an increase times a smaller decrease is a decrease, and . = since the identity may be obtained by the multiplication of an increase and a decrease one the exact inverse of the other.

δ δ δ • ∆3 δ < < {<} < > {<, . =, >} < . = {<, . =} > < {<, . =, >} > > {>} > . = { .

=, >}

The table for ma k ∆6 = • ∆6 can be computed analogously. Note however that there exist partitions ∆ such that ma k ∆ = • ∆ . For instance, we can choose ∆ = {I, ¬I} where I = {(1, 1)} and ¬I = R 2 + \ I. In this case we have for any k = 1 that ma k ∆ (I, I) = {(k, k)} ∆ = {¬I} = I • ∆ I, while for k = 1 it holds that ma 1 ∆ (I, I) = {I} = I • ∆ I. This partition ∆ also illustrates that difference abstractions ma k ∆ may depend on the choice of parameter k.

Example 2. We consider Michaelis-Menten laws that are enhanced with an additional activator, mm k1,k2 (subs: d, act:

d ) = k1 • d • d k2+d
where k, k ≥ 0. It can be shown that (mm k1,k2 ) ∆3 = • ∆3 and (mm k1,k2 ) ∆6 = • ∆6 . Therefore, these difference abstractions are again independent of the choice of the parameters. Furthermore, for all parameters k, k 1 , k 2 > 0, the two abstractions with partitions ∆ ∈ {∆ 3 , ∆ 6 } satisfy:

ma k ∼ ∆ mm k1,k2
This shows that these abstractions cannot distinguish between the mass-action laws and the enhanced Michaelis-Menten laws. It should be noticed though that there exist difference abstractions with other ∆'s that are able to do so.

It should be noticed that R ∆ is always a finite relation, since ∆ is chosen to be finite. The relation R in contrast, may contain infinitely many tuples. As a consequence, infinitely much information may be abstracted away, in particular the choice of the parameters of the kinetic functions.

The information preserved when abstracting with respect to ∆ 3 is able to distinguish between inhibitors and activators, since these correspond to decreases > and respectively increases <, while the distinction between activators and accelerators is very faint. In contrast, the refined difference abstraction with respect to ∆ 6 is able to distinguish them clearly. Activation corresponds to increases from zero ⇑ and acceleration to increases in ↑ ∪ ⇑. Similarly, the refined abstraction can clearly distinguish reaction knockouts from inhibitions that only slow down the reaction. The former corresponds to knockout decreases ⇓ of an active reaction and the latter to slow-down decrease ↓ of an active reaction. Finally, removals of reactions from the network correspond to difference relation ≈, while ∼ corresponds to the no-change of an active reaction.

Part II

Prediction Methods

We now present prediction methods for network changes, that apply to all models defined in our language.

We show in Section 6 how to derive difference constraints from steady state equations based on abstract interpretation. In Section 7 we explain the main ideas of how difference constraints can be used for qualitative reasoning about reaction networks, as needed for our prediction algorithms. In Section 8 we show how to infer from the solution set of a difference constraint the change predictions that it implies, and whether these predictions are just compatible for the target or even safe. Our tools for modeling reaction networks with partial kinetic information and prediction of network changes is presented in Section 9.

Abstract Interpretation to Difference Constraints

We now explain how to interpret steady state equations abstractly as difference constraints, which are used to apply qualitative reasoning on reaction networks modeled in our language, as we will show in the next section.

The idea is to lift the difference abstraction . ∆ from relations over R + to relations over ∆ to the level of constraints defining such relations. For instance, the arithmetic equation x A = ma k (z A , z B ) can be abstracted to a difference constraint that defines the relation ma k ∆ . We write this difference constraint as x A ∈ ma k (z A , z B ), since now the variables are interpreted by values of ∆ and ma k is interpreted as the set valued function ma k ∆ . It should be noticed that the relation ma k ∆ is finite and independent of the unknown parameter k, i.e. the unknown parameter has been abstracted away successfully.

Arithmetic Constraints

Arithmetic constraints were used to define the steady state semantics of reaction networks. More formally, an arithmetic constraint is a conjunctive logic formula with existential quantifiers with the following abstract syntax:

φ ::= x=κ (i) (x 1 , . . . , x k ) | x = x 1 + x 2 | x = x 1 • x 2 | x 1 = x 2 | x ∈ S | φ ∧ φ | ∃x.φ
where i ∈ N, κ : R k + → R + , x, x 1 , x 2 are variables, and S ⊆ R + a finite set. The expression κ (i) is a variable for a kinetic function that is similar to κ, i.e. an implicitly existentially quantified variable that is subject to the similarity constraint κ (i) ∼ κ.

A solution of an arithmetic constraint φ with n variables can be identified with a tuple in R n + since we assumed a total order on the variables. Therefore, the solution set sol (φ) of a formula φ satisfies sol (φ) ⊆ R n + .

Lemma 1. The steady state equations of any reaction network N can be rewritten in linear time to an equivalent arithmetic constraint φ N so that:

sol (φ N ) = E N Proof.
It is sufficient to flatten the terms in steady state equations, by introducing new existentially bound variables for all nested subterms. The resulting constraint φ N is clearly equivalent, and thus it satisfies sol (φ N ) = E N by Definition 2 of the exchange relation.

Difference Constraints

A difference constraint is a conjunctive logic formula with existential quantifiers with the following abstract syntax, where x is a variable, d ∈ R + , S ⊆ R + , and κ : R k

+ → R + : difference relation t ::= x | d set of difference relations s ::= t | S | s + s | s • s | κ(s 1 , . . . , s k ) difference constraints ψ ::= t ∈ s | t=t | ψ ∧ ψ | ∃x.ψ
The intuition of the semantics of a difference constraint is that its variables denote difference relations of some set ∆, while all real-valued relations in the difference constraint are ∆-abstracted to + ∆ , • ∆ , S ∆ , and κ ∆ respectively. A real number d ∈ R + is interpreted as the unique element of {d} ∆ . More formally, the semantics of a difference constraint for a partition ∆ of difference relations is defined in Figure 10. It assumes a variable assignment α into ∆. The semantics of a term t is the difference relation t ∆,α ele ∈ ∆, the semantics of a term s is the set of difference relations s ∆,α set ∈ 2 ∆ , and the semantics of a difference constraint ψ is the Boolean ψ ∆,α ∈ B. Since any term of type t is also term of type s, it has two different semantics. These satisfy t ∆,α set = { t ∆,α ele }. Whether a term of t is given the element or set semantics depends on the context in which is used.

Difference relations:

x ∆,α ele = α(x) d ∆,α ele = unique element of {d} ∆ Sets of difference relations: The set of all ∆-solutions of a difference constraints ψ with free variables x 1 , . . . , x n -listed in their order -is the following subset of ∆ n :

t ∆,α set = { t ∆,α ele } S ∆,α set = S ∆ s • s ∆,α set = s ∆,α set • ∆ s ∆,α set s + s ∆,α set = s ∆,α set + ∆ s ∆,α set κ(s 1 , . . . , s k ) ∆,α set = κ ∆ ( s 1 ∆,α set , . . . , s k ∆,α set ) Difference constraints: t ∈ s ∆,α = ( t ∆,α ele ∈ s ∆,α set ) t = t ∆,α = ( t ∆,α ele = t ∆,α ele ) ψ ∧ ψ ∆,α = ψ ∆,α ∧ B ψ ∆,α ∃x.ψ ∆,α = (∃δ ∈ ∆. ψ ∆,α[x/δ] = 1)
sol ∆ (ψ) = {(α(x 1 ), . . . , α(x n )) | ψ ∆,α = 1}

Abstract Interpretation

We can now abstract from arithmetic constraints by interpreting them as difference constraints:

x=κ (i) (x 1 , . . . , x k ) = df x ∈ κ(x 1 , . . . , x k ) x ∈ S = df x ∈ S x = x 1 + x 2 = df x ∈ x 1 + x 2 φ ∧ φ = df φ ∧ φ x = x 1 • x 2 = df x ∈ x 1 • x 2 ∃x.φ = df ∃x. φ x 1 =x 2 = df (x 1 =x 2 )
An important point here is that the variables κ (i) for the partially known kinetic functions are replaced by the the kinetic functions κ (which are interpreted as the κ ∆ in the semantics). Since the kinetic functions such as exp and ma have definitions by arithmetic expressions containing calls of the basic functions Inh and Acc and some constants only, their applications can be replaced by arithmetic expressions (interpreted over ∆), in which the only remaining kinetic functions will be Inh, Acc, and the constants d ∈ R + . In this way, the simplified steady state equations in Figure 6 for the basic P Ilv-Leu network are abstracted to the difference constraints in Figure 11.

For instance, the abstraction of x = ma (i) (subs:x 1 ) becomes the equation x = x 1 , and the abstraction of x = ma (i) (subs:x 1 , acc: It should also be noticed that a membership constraint, such as x ∈ {0, 1}, remains unchanged syntactically by abstract interpretation, but afterwards it is interpreted semantically over ∆, so that it becomes for example equivalent to x ∈ {0, 1}

x 2 ) the membership constraint x ∈ x 1 • Acc(x 2 ). For x = exp (i) (subs: x 1 , inh: x 2 , inh:x 3 ) the abstraction yields x ∈ x 1 •Inh(x 2 +x 3 ). The equation x = exp (i) () is interpreted abstractly to x ∈ 1. Note that 1 ∆3,α set = { . =} while 1 ∆6,α set = {∼} for any variable assignment α. v r2 ∈ v r 2 + y Leu y Leu = z Leu v r1 ∈ z CodY • Inh(z CcpA ) v r1 = z BS CodY v r2 ∈ z BSCcpA • Inh(z BS CodY + z Leu +z BSTnrA ) v r 2 = z P Ilv-Leu v r7 = z TnrA v r7 = z BSTnrA y Leu = z P Ilv-Leu v r9 ∈ z CcpA • Inh(z BSTnrA ) v r9 = z BSCcpA v r14 ∈ 1 v r14 = z CcpA x CodY = z CodY x TnrA = z TnrA
∆ . Over ∆ 6 this is equivalent to x ∈ {⇑, ⇓, ∼, ≈}, while over ∆ 3 it is always true.

Theorem 1 (Soundness). For any arithmetic constraint φ and partition ∆ ⊆ 2 R+×R+ into difference difference relations, the set of ∆-solutions of the abstract interpretation φ over-approximates the ∆-difference abstraction of the set of solutions of φ, that is:

sol (φ) ∆ ⊆ sol ∆ ( φ ).
Proof. We first note that for any two relations

R 1 , R 2 ⊆ R n + : (intersect) (R 1 ∩ R 2 ) ∆ ⊆ R 1 ∆ ∩ R 2 ∆
This can be verified straightforwardly from the definition of the difference abstraction. It will turn out to be the reason why proper approximations may appear when abstracting conjunctions.

Then we proceed by induction on the structure of arithmetic constraints φ. For the variables in the proof, we always assume for simplicity that they are ordered x 1 x 2 . . . x k x.

Case φ is x = κ (i) (x 1 , . . . , x k ). By definition of solutions we have that sol (φ) = {(d 1 , . . . , d k , d) ∈ R k+1 + | d = κ (d 1 , . . . , d k ), κ ∼ ∆ κ}, i.e., sol (φ) = {κ | κ ∼ ∆ κ}. Thus, sol (φ) ∆ = κ ∆ . Furthermore, φ is equal to x∈κ(x 1 , . . . , x k ) so that sol ∆ ( φ ) = κ ∆ = sol (φ) ∆ . Case φ is x = x 1 + x 2 . This follows, since sol (φ) ∆ = + ∆ . Furthermore φ is equal to x ∈ x 1 +x 2 so that sol ∆ ( φ ) = + ∆ and thus sol (φ) ∆ = sol ∆ ( φ ). Case φ is x = x 1 • x 2 .
This follows in analogy to the case of addition.

Case φ is

x 1 = x 2 . We first note that {(d, d) | d ∈ R + } ∆ = {(δ, δ) | δ ∈ ∆} since ∆ is a partition of R + × R + . We now can conclude as follows: sol (φ) ∆ = {(d, d) | d ∈ R + } ∆ = {(δ, δ) | δ ∈ ∆} = sol ∆ φ .
Case φ is x ∈ S. In this case we have sol (φ) = S, so that sol (φ)

∆ = S ∆ = sol ∆ (x ∈ S) = sol ∆ ( φ ).
Case φ is φ 1 ∧ φ 2 . Hence sol (φ) = sol (φ 1 )∩sol (φ 2 ) and thus as claimed by (intersect) at the beginning sol (φ)

∆ ⊆ sol (φ 1 ) ∆ ∩ sol (φ 2 ) ∆ .
From the induction hypothesis applied to φ 1 and φ 2 , it follows that sol (φ i ) ∆ ⊆ sol ( φ i )

for i = 1, 2, and thus sol (φ)

∆ ⊆ sol (φ 1 ) ∆ ∩ sol (φ 2 ) ∆ ⊆ sol ∆ ( φ 1 ) ∩ sol ∆ ( φ 2 ) = sol ∆ ( φ ).
Case φ is ∃x.φ . There case where x does not occur freely in φ is easy, since we can drop the quantifier ∃x. Otherwise, we can assume that φ has the free variables x 1 , ..., x k and that x = x i where 1 ≤ i ≤ k. For any set D and tuple set S ⊆ D k we define the projection

Π x (S) = {(d 1 , . . . , d i-1 , d i+1 , . . . , d k ) | (d 1 , . . . , d k ) ∈ S}.
We first note that:

(project) Π x (sol (φ )) ∆ = (Π x (sol (φ ))) ∆ All elements of R + × R + belong to some difference in ∆.
Based on the equation on projection, we can conclude as follows:

sol (∃x.φ ) ∆ = (Π x (sol (φ ))) ∆ = project Π x (sol (φ )) ∆ ⊆ ind.hyp. Π x sol ( φ ) = sol (∃x. φ ) = sol ∆ ( ∃x.φ ).
This theorem shows for any reaction network N with steady state equations equivalent to φ N that the set of ∆-solutions of the abstract interpretation φ N is a correct over-approximation of the abstraction E N ∆ of the exchange relation of N .

Corollary 1. If E N = sol (φ N ) then (E N ) ∆ ⊆ sol ∆ ( φ N ).
Proof. This follows immediately from Theorem 1.

Qualitative Reasoning with Difference Constraints

We next illustrate how to reason qualitatively about reaction networks by investigating solution sets of difference constraints. Since difference constraints have finite domains, we can compute all solutions of difference constraints by using finite domain constraint programming. In simple cases, we can use constraint simplification for this purpose, similarly to what is done by the propagation mechanisms for finite domain constraint solvers.

Specifying Change Targets

The first thing we have to do is to specify a target for a network changes. For instance, me may to change the P Ilv-Leu network in order to increase the outflux of Leu. In this case, the overproduction target can be expressed by the formula y Leu = <. It should be notice that this formula is not a difference constraint. Indeed, it can only be interpreted over ∆ 3 but not over arbitrary ∆'s. In general, a change target is some ∆-difference formula for some fixed domain ∆:

= . = . = no 7. > < . = no 8. < > . = no 9. . = < > no 10. < . = > no 11. > < > no 12. < > > no 13. < < > no Figure 12:
The set of all solutions of the difference constraint in Figure 11, derived from the basic PIlv-Leu network (without OP BkL-Bcd ), and their interpretation with respect to the target y Leu = <.

Definition 4. The set of ∆-difference formulas Ψ is the least set that contains all difference constraints φ, all formulas x ∈ ∆ where x is a variable and ∆ ⊆ ∆, and that is closed under the first-order connectives, ie.:

Ψ ::= φ | x ∈ ∆ | Ψ ∧ Ψ | ¬Ψ | ∃x.Ψ
We define the ∆-difference formula x = δ as syntactic sugar for x ∈ {δ} for any δ ∈ ∆ and variable x.

Constraint Solving

We illustrate how to use constraint solving to answer the question from Section 4.3: which changes of the influxes of the basic P Ilv-Leu network may cause the change target y Leu = <? In order to find the answer, we interpret the difference constraint inferred for the network in Figure 11 over ∆ 3 , which is good enough for this simple case.

We now consider the difference constraint in Figure 12, which is derived from the basic PIlv-Leu network (without OP BkL-Bcd ). The control variables of this network are x CodY and x TnrA , since these represent the changes that we want to control. Besides the control variables, we are also interested in the free variables of the target formula, which in our case correspond only to y Leu . All the other variables of the constraint are considered as local, so that their values are projected away in solution sets.

The set of all solutions of this constraint could be enumerated by our finite domain constraint solver. There are 13 solutions after projection to the global variables x CodY , x TnrA and y Leu . Only 5 solutions do also satisfy the target y Leu = <. From those, only the first (inh+) Inh(s) • Inh(s ) ⇒ Inh(s + s ) three are safe with respect to the target, in that no other solution with the same values of x CodY and x TnrA violates the target. This shows that the target will be safely satisfied if we either decrease the influx of CodY, the influx of TnrA, or both of them.

(no 1 ) Inh(1) ⇒ 1 (no 3 ) t • 1 ⇒ t (ip) x + x ⇒ x (no 2 ) Acc(1) ⇒ 1 (no 4 ) 1 • t ⇒ t (si) t ∈ t ⇒ t = t (bv) ∃x. (x = t ∧ ψ) ⇒ ψ[t/x] (inh) t ∈ Inh(t + s) ⇒ t ∈ Inh(s)

Constraint Simplification

Since the P Ilv-Leu network is quite simple, one can obtain the same predictions based on constraint simplification by term rewriting without any search. This kind of formal reasoning is very similar to the intuitive reasoning from Section 4.3.

Here we present a simplification algorithm that is correct only for specific choices of ∆, as for instance ∆ 3 and ∆ 6 . The first restriction we need is that Acc ∆ and Inh ∆ are functional relations, so that we can rewrite x ∈ Inh(t) to x = Inh(t) and then eliminate x by substitution with Inh(t).

In order to do so, we consider a variant of difference constraints in which applications Acc(t) and Inh(t) are used for defining a single difference relation, rather than a set thereof. difference relation

t ::= x | d | Acc(t) | Inh(t) set of difference relations s ::= t | S | s + s | s • s | Acc(s) | Inh(s) difference constraints ψ ::= t ∈ s | t=t | ψ ∧ ψ | ∃x.ψ
In Figure 13 we present a collection of basic term rewrite rules for simplifying difference constraints. Their correctness under the above assumptions on ∆ is straightforward. Rule (bv) replaces equal by equal while eliminating existentially bound variables (all variables z A and v ri are implicitly existentially quantified). The simplification rules (no i ) remove the constant 1 (which over ∆ 3 is interpreted as .

= and over ∆ 6 as ∼). The third rule (si) simplifies membership in singletons to equality. Rule (ip) expresses the idempotence of addition. The result of the basic simplification for the difference constraints for the basic P Ilv-Leu network with the rewrite rules in Figure 13 is the following constraint:

y Leu ∈ Inh(x TnrA ) • Inh(x CodY + y Leu + x TnrA ) .
The constraint can be further simplified by applying the rewrite rule (inh+) in Figure 14, whose correctness requires a further restriction on ∆:

y Leu ∈ Inh(x TnrA + x CodY + y Leu + x TnrA ) .
With a further application of simplification rule (ip), we obtain:

y Leu ∈ Inh(x CodY + y Leu + x TnrA ) .
When working over ∆ 3 and assuming our target y Leu = < then we can simplify the constraint further to: < ∈ Inh(x CodY + x TnrA ) This is equivalent to

x CodY = > ∨ x TnrA = > .
This can be satisfied by decreasing the influx of either CodY or TnrA. Thus, we obtain the same safe solutions as before.

Predicting Multiple Changes

Some changes of influxes or reaction knockouts may have more than one possible outcome, so that the target get either get satisfied or not -while others may have only one outcome. This observation was made already in Section 7 with the solutions that were safe or unsafe for leucine overproduction.

For example, reconsider Figure 12 with the solutions of the difference constraint for the basic P Ilv-Leu network. Consider the solution 5, 8, and 12. All three solutions require an increase of the CodY influx and a decrease of the TnrA influx at the same time. But the outcomes are different: solution 5 satisfies the target of leucine overproduction, while solutions 8 and 12 do not since they decrease the outflux of Leu or respectively leave it unchanged. Therefore, we said that solution 5 is not safe for the target. On the other hand, some changes may have only one possible outcome, as for safe solutions 1, 2 and 3.

It should be noticed that multiple outcomes are possible because of the lack of information on kinetic parameters. In the example, this does not let us determine whether the effect of the increase of CodY influx prevails on the decrease of TnrA, or the other way around, or whether they cancel each other out.

When looking for the overproduction of metabolites it is then important to distinguish the safe changes that for sure produce the desired effect because the overproduction happens for all the possible outcomes, from the compatible changes that simply may produce the target effect, because for some other outcomes the effect does not happen.

Definition 5 (N-fold Changes). The set of control variables of a network N = (S, R, I, O, K) is {x s | s ∈ I} ∪ {o r | r ∈ K}. A ∆-control assertion for N is a ∆-difference formula of the following form:

x 1 = δ 1 ∧ x 2 = δ 2 ∧ • • • ∧ x p = δ p
where x 1 , . . . , x p are pairwise distinct control variables of N and δ 1 , . . . , δ p ∈ ∆. A control assertion is called total if it contains all control variables of N and partial otherwise. A n-fold change is a total control assertion in which the number of variables with values in ∆ \ { d ∆,α ele | d ∈ R + } is equal to n. For ∆ 6 , an inflow change has the form x s = ↑ or x s = ↓ where s ∈ I, and a knockout has the form o r = ⇓ for some knockout candidate r ∈ K. When inflow changes and knockouts are combined, we talk about mixed changes, as opposed to strict changes, when only knockouts or only inflow changes are involved. Several changes are generally possible for a given network, involving one or more inflow variables as for the basic P Ilv-Leu network, but possibly also one or more knockout variables as discussed later for the refined P Ilv-Leu network.

Definition 6 (Safeness). Let N be a network and ∆ a set of difference relations, Θ be a be a ∆-difference formula, and Ψ be a ∆-change. The change Ψ is said safe for a target Θ if φ N ∧ Ψ |= Θ (over ∆). We say that Ψ is unsafe for

Θ if φ N ∧ Ψ ∧ Θ is satisfiable but φ N ∧ Ψ |= Θ.
So for example there are two single ∆ 3 -changes and one double ∆ 3 -change in the basic P Ilv-Leu network, that are safe with respect to the target of leucine overproduction. Note that this network has no knockout candidates, so all changes there are inflow changes. The two single changes correspond to solutions 1 and 2 in Figure 12 and the double change to solution 3. Remarkably, this double safe change is obtained as a combination of two single safe changes. Conversely, when a safe change is combined with an unsafe change or with a change which is not compatible with the target, the result is often unsafe as well, as for solution 4 obtained by the combination of 1 and 9. According to our experience, the combination of safe or unsafe changes with other unsafe changes or changes incompatible with the target gives usually rise to unsafe changes. Exceptions are however possible as discussed in Section 8.1.

Depending on the complexity of the network there may exist no safe change at all.

Application example

The refined P Ilv-Leu network in Figure 7 allows several changes, given by all the combinations of inflow changes discussed for the simple P Ilv-Leu network together with the knockouts of r 1 , r 3 , r 7 , r 14 . All the safe and unsafe changes are listed respectively in Appendix 13.1 and 13.2. The presence of the species OP BkL-Bcd influenced by CodY and TnrA modifies considerably the effect of the changes discussed in the simple network. There the two inflows acted as simple inhibitors of Leu, while in the refined network both the effects of inhibition and acceleration are present and make unpredictable the outcome of any strict inflow change.

In fact, all the strict inflow changes are compatible with leucine increase but unsafe, as shown by solutions 1-4 and 20-23 in Appendix 13.2. The simplest safe changes are instead represented by the knockouts of r 1 , r 3 or r 7 , corresponding to solutions 1-3 in Appendix 13.1. These can be combined at will to obtain more complex changes, that are still safe as shown by solutions 4-6 and 11. Intuitively, each of these knockouts blocks one of the three inhibitions acting on the production of leucine, either indirectly through r 2 or directly on r 8 . Conversely, the knockout of r 14 alone is not compatible with leucine increase, as one would expect since the only effect of this reaction is an acceleration of the production of the unique Leu precursor.

Further safe changes can be remarkably obtained by combining inflow changes (that are unsafe when applied alone, as previously discussed, or even incompatible with leucine outflow increase) with targeted knockouts. For example, the knockout of r 1 can block the inhibitory effect of CodY on Leu, so that an increase of CodY inflow also contributes to the overproduction of leucine. Intuitively, in Figure 7 we can see that once r 1 is blocked, an increase of CodY inflow has the only effect of inhibiting the inhibition of the production of Leu carried out by OP BkL-Bcd , whose net effect is an acceleration of leucine outflow. In the same way r 7 can be knocked out and the inflow of TnrA can be increased to obtain the safe mixed double change 10 in Appendix 13.1.

Conversely, it is possible to block the positive effect of CodY and TnrA on Leu by knocking out r 3 , getting back in practice the basic P Ilv-Leu network. Indeed, if r 3 is knocked out we obtain solutions analogous to those listed in Figure 12, corresponding respectively to safe changes 8, 9, 20 in Appendix 13.1 and unsafe changes 47, 48 in Appendix 13.2. The presence of the knockout of r 3 allows however several other unsafe changes obtained by combination with CodY and TnrA inflow changes.

Safe double mixed changes can be further combined to obtain more complex safe changes, represented by solutions 12-19, 21-33 in Appendix 13.1.

Application to more complex networks

The basic and refined P Ilv-Leu networks are part of a more complex network module presented in [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF], where single mixed knockouts were discussed. With respect to leucine overproduction, this complex network is an example of a model for which there exists no safe single or double knockout. Despite the complexity of the network, qualitative reasoning allowed us to exclude 7 over 21 knockout candidates because they turned out to be incompatible with leucine overproduction. This shows how qualitative reasoning can be useful even in the absence of safe solutions, when uncertainty allows reasoning only about compatibility or incompatibility of changes with respect to the desired target.

Work on improving the predictions for the complex network in [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] are in progress. Preliminary results show us the presence of safe changes appearing after a minimum of three knockouts, remarkably in agreement with the choice of the 6 knockouts previously tested in wet lab experimentation.

Modeling and Prediction Tool

We have implemented a tool for modeling reaction networks in our language and performing change predictions. The input of the tool are a reaction network in XML syntax and an overproduction target. From that, the steady state semantics is derived, from which the corresponding difference constraints are computed and solved. The set of all solutions is then analyzed in order to filter out solutions that are compatible with the target or even safe. XML Syntax. As its input, our tool uses an Xml syntax for reaction networks that is then converted into an equivalent graphical syntax by an Xslt transformation, which produces latex graphics using the Tikz package. These graphics are included into the sources of the present article. For instance, the Xml syntax for the basic P Ilv-Leu network from Figure 4 is given in Appendix 14.

Since both syntaxes are equivalent, the Xml syntax could be generated from the graphical syntax. There exist tools for graphical input for Petri nets, but their sources are not freely available. And as it turned out for our purposes, the XML syntax was by far good enough as input syntax.

Constraint Solver. Difference constraint can be solved over a fixed domain ∆ by finite domain constraint solving. The previous constraint solver from [START_REF] John | Knockout Prediction for Reaction Networks with Partial Kinetic Information[END_REF] was implemented in Scala [START_REF] Odersky | The evolution of scala: Ple'14 keynote[END_REF] from scratch. This solver could enumerate the n-best solutions where n ≤ 3000, while consuming considerable time with more than 10 minutes.

We now implemented a new solver by compiling difference constraints over ∆ to MiniZinc [START_REF] Nethercote | Minizinc: Towards a standard CP modelling language[END_REF]. Thereby the difference constraints can be solved much more efficiently. Therefore it can indeed return the complete set of all solutions sets in all our application, while enumerating more than 5000 solutions in less then a second. An example for the MiniZinc constraint formulation of the difference constraint in Figure 11 over ∆ 6 is given in the Appendix 15.

Analysis of Solutions Sets. A considerable effort was also spent for a proper analysis of solutions sets. Most importantly, they must be presented in a readable output format, since they tend to be very large. For this, we sorted and grouped solution sets, eliminated duplicates and removed irrelevant variables. We also compute the sets of safe solutions automatically, based on XPath queries with data joins. All these transformations were also done in Xslt.

Conclusions

We have presented a formal modeling language for chemical reaction networks with partial kinetic information, and shown how to abstract away from the unknowns thanks abstract interpretation. We have illustrated that this allows us to reason qualitatively about such networks, so that we can predict influx changes, reaction knockouts, or multiple change of these types. We have illustrated the relevance of the notion of safe changes, and how to compute them by inspecting the set of all solutions of a difference constraint.

An immediate question for future work is to find safe multiple changes for big reaction networks as the one in [START_REF] Coutte | Modeling Leucine's Metabolic Pathway and Knockout Prediction Improving the Production of Surfactin, a Biosurfactant from Bacillus Subtilis[END_REF] and to validate them in the wet lab. An important question for future work on the longer scale is how to develop finer abstractions for quantitative predictions. <context id="100" x="21.5" y="-1.75"> <input spec="CodY"/> </context> <context id="101" x="21.5" y="-6.25"> <input spec="TnrA"/> </context> <context id="102" x="11.05" y="-6.25"> <output spec="Leu"/> </context> <!--reactions --> <reaction id="1" x="19.00" y="-2.75" comment="bind \M\CodY to \M\PIlvLeu for inhibition"> <kinetics id="exp" angle="-120"/> <inhibitor spec="CcpA"/> <activator spec="CodY"/> <product spec="BSCodY"/> </reaction> <reaction id="1'" x="15.0" y="-2.05" comment="implicit degradation"> <kinetics id="ma" angle="45"/> <reactant spec="BSCodY"/> </reaction> <reaction id="2" x="14.20" y="-4.00" comment="activate \M\PIlvLeu promoter"> <kinetics id="exp" angle="30"/> <inhibitor edgecluster="BSCodYa"/> <accelerator spec="BSCcpA"/> <inhibitor spec="Leu"/> <inhibitor edgecluster="BSTnrAa"/> <product spec="PIlv-Leu"/> </reaction> <reaction id="2'" x="11.8" y="-3.3" comment="implicit degradation"> <kinetics id="ma" angle="30"/> <reactant spec="PIlv-Leu"/> </reaction> <reaction id="3" x="19.0" y="-7.25" comment="bind \M\BkdR to \M\BkL \M\Bcd promoter"> <kinetics id="exp" angle="-30"/> <inhibitor spec="TnrA"/> <inhibitor edgecluster="CodY2b"/> <product spec="OpBkLBcd"/> </reaction> <reaction id="3'" x="13.65" y="-7.25" comment="implicit degradation"> <kinetics id="ma" angle="30"/> <reactant spec="OpBkLBcd"/> </reaction> <reaction id="7" x="19.00" y="-5.25" comment="bind \M\TnrA to \M\PIlvLeu promoter for inhibition"> <kinetics id="exp" angle="120"/> <activator spec="TnrA"/> <product spec="BSTnrA"/> </reaction> <reaction id="7'" x="15.0" y=" -4.7" comment="implicit degradation"> <kinetics id="ma" angle="30"/> <reactant spec="BSTnrA"/> </reaction> <reaction id="8" x="12.5" y="-5.00" comment="\M\PIlvLeu expression followed by \M\Leu production"> <kinetics angle="45" id="exp"/> <inhibitor spec="OpBkLBcd"/> <reactant spec="PIlv-Leu"/> <product spec="Leu"/> </reaction> <reaction id="9" x="17.5" y=" -4.0" comment="bind \M\CcpA to \M\PIlvLeu promoter without \ M\BSTnrA loop"> <kinetics id="exp" angle="-30"/> <activator spec="CcpA"/> <inhibitor spec="BSTnrA"/> <product spec="BSCcpA"/> </reaction> <reaction id="9'" x="15.0" y="-3.3" comment="implicit degradation"> <kinetics id="ma" angle="30"/> <reactant spec="BSCcpA"/> </reaction> <reaction id="14" x="20.5" y="-4.0" comment="expressions of CcpA"> <kinetics id="exp" angle="30"/> <product spec="CcpA"/> </reaction> <reaction id="14'" x="17.5" y="-3.3" comment="implicit degradation"> <kinetics id="ma" angle="30"/> <reactant spec="CcpA"/> </reaction> <reaction id="15'" x="17.5" y="-1.75" comment="\CodY deactivation"> <kinetics id="ma" angle="30"/> <reactant spec="CodY"/> </reaction> <reaction id="16'" x="17.5" y="-6.25" comment="\TnrA deactivation"> <kinetics id="ma" angle="30"/> <reactant spec="TnrA"/> </reaction> </network>

o r1 o r3 o r7 1. ⇓ ∼ ∼ 2. ∼ ⇓ ∼ 3. ∼ ∼ ⇓ o r1 o r3 o r7 4. ⇓ ⇓ ∼ 5. ⇓ ∼ ⇓ 6. ∼ ⇓ ⇓ o r1 o r3 o r7 x CodY x TnrA 7. ⇓ ∼ ∼ ↑ ∼ 8. ∼ ⇓ ∼ ↓ ∼ 9. ∼ ⇓ ∼ ∼ ↓ 10. ∼ ∼ ⇓ ∼ ↑ o r1 o r3 o r7 11. ⇓ ⇓ ⇓ o r1 o r3 o r7 x CodY x TnrA 12. ⇓ ⇓ ∼ ↓ ∼ 13. ⇓ ⇓ ∼ ∼ ↓ 14. ⇓ ⇓ ∼ ↑ ∼ 15. ⇓ ∼ ⇓ ∼ ↑ 16. ⇓ ∼ ⇓ ↑ ∼ 17. ∼ ⇓ ⇓ ↓ ∼ 18. ∼ ⇓ ⇓ ∼ ↓ 19. ∼ ⇓ ⇓ ∼ ↑ o r3 x CodY x TnrA 20. ⇓ ↓ ↓ o r1 o r3 o r7 x CodY x TnrA 21. ⇓ ⇓ ⇓ ↓ ∼ 22. ⇓ ⇓ ⇓ ∼ ↓ 23. ⇓ ⇓ ⇓ ∼ ↑ 24. ⇓ ⇓ ⇓ ↑ ∼ o r1 o r3 o r7 x CodY x TnrA 25. ⇓ ⇓ ∼ ↓ ↓ 26. ⇓ ⇓ ∼ ↑ ↓ 27. ⇓ ∼ ⇓ ↑ ↑ 28. ∼ ⇓ ⇓ ↓ ↓ 29. ∼ ⇓ ⇓ ↓ ↑ o r1 o

Figure 1 :

 1 Figure 1: A reaction network describes a system of chemical reactions that acts on a chemical solution, the inflow of molecules from the context, and the outflow of molecules into the context.

Figure 2 :

 2 Figure 2: Steady state equations of a reaction network N = (S, R, I, O, K).

Definition 2 .

 2 Any reaction network N = (S, R, I, O, K) with n = #I inflow species, m = #dom(O) outflow species, and o = #K knockout candidates defines an exchange relation E N ⊆ R n+m+o +

Figure 3 :

 3 Figure 3: An enriched reaction with a partially known rate law ∼κ . It has substrate S twice, inhibitor I, accelerator A , activator A, and product P once.

  (), which means that the kinetic function of r 14

Figure 4 :

 4 Figure 4: Basic reaction network of the regulation of promoter PIlv-Leu in B. subtilis without knockout candidates.

Figure 5 :

 5 Figure 5: Steady state equations for the basic reaction network in Figure 4.

Figure 6 :

 6 Figure 6: Simplified steady state equations for the basic PIlv-Leu network.

Figure 7 :

 7 Figure 7: Refined PIlv-Leu network with regulation of Leu production by the BkL-Bcd operon. As knockout candidates, we added the five reactions in light orange.

  2a. decrease the influx of CodY and knockout reaction r 3 , or 2b. increase the influx of CodY and knockout reaction r 1 , or 2c. decrease the influx of TnrA and knockout reaction r 3 , or 2d. increase the influx of TnrA and knockout reaction r 1 .

Figure 8 :

 8 Figure 8: The difference relations of partition ∆ 3 = {<, >, . =}.

Figure 9 :

 9 Figure 9: The difference relations of partition ∆ 6 = {↑, ⇑, ↓, ⇓, ∼, ≈}.

Definition 3 .

 3 The general abstraction applies to relations R ⊆ R p + . This subsumes the case of kinetic function κ of arity p -1 and of exchange relations R N for networks N with p = n + m + o. We define ∆-difference abstraction of R as follows:R ∆ = {(δ 1 , . . . , δ p ) ∈ ∆ p | exist (d 1 , d 1 ) ∈ δ 1 , . . . ,(d p , d p ) ∈ δ p : (d 1 , ..., d p ) ∈ R and (d 1 , ..., d p ) ∈ R} So for instance, consider the exchange relation E N for some reaction network N . Its difference abstraction E N ∆ then expresses how the tuples in E N may change when moving from one steady state of N to another. Two kinetic functions κ 1 , κ 2

Figure 10 :

 10 Figure 10: Semantics of difference constraints over ∆.

Figure 11 :

 11 Figure 11: Difference constraints for the basic PIlv-Leu network (without OP BkL-Bcd ), inferred from the simplified steady state equations in Figure 6.

Figure 13 :

 13 Figure 13: Basic simplification rules.

Figure 14 :

 14 Figure 14: Special simplification rule

Part III Appendix 13 .

 13 Solutions for Leucine Increase 13.1. Safe Changes for ∆ 6

x

  CodY x TnrA y Leu

					target	quality
					satisfied?
	1. 2.	> . =	. = >	< <	yes yes	safe safe
	3.	>	>	<	yes	safe
	4.	>	<	<	yes	unsafe by 7. or 11.
	5.	<	>	<	yes	unsafe by 8. or 12.
	6.

.

  r3 o r7 x CodY x TnrA

	o r1 o r3 o r7 o r14 o r1 o r3 o r7 o r14 x CodY x TnrA o r1 o r3 o r7 o r14 x CodY x TnrA
	58. ⇓ 89. ⇓	24. ⇓ ⇓ ⇓	∼ ∼		⇓	⇓ ⇓	∼	⇓ ↓ ↓	∼ ↓
	59. ⇓ 90. ⇓	25. ⇓ ⇓ ⇓	∼ ∼		∼	⇓ ⇓	⇓	⇓ ∼ ↓	↓ ↑
	60. ⇓ 91. ⇓	26. ∼ ⇓ ⇓	∼ ∼		⇓	⇓ ⇓	⇓	⇓ ∼ ↑	↑ ↓
	61. ⇓ 92. ⇓ o r1 o r3 o r7 o r14 x CodY x TnrA ⇓ ∼ ⇓ ↑ ∼ ⇓ ∼ ⇓ ↑ ↑ 62. ⇓ ∼ ⇓ ⇓ ↓ ∼ 93. ⇓ ∼ ⇓ ⇓ ↓ ↓
	27. ⇓ 63. ⇓ 94. ⇓		⇓ ∼ ∼			∼ ⇓ ⇓			∼ ⇓ ⇓		∼ ∼ ↓	↑ ↓ ↑
	30. ⇓ 28. ⇓ 64. ⇓ 95. ⇓ 31. ⇓ 29. ⇓ 65. ⇓ 96. ⇓ 32. ⇓ 30. ⇓ 66. ∼ 97. ∼ 33. ⇓ 32. ⇓ 31. ⇓ 67. ∼ 98. ∼	⇓ ⇓ ⇓ ∼ ∼ ∼ ⇓ ∼ ⇓ ∼ ⇓ ∼ ⇓ ∼ ⇓ ∼ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ∼ ∼ ∼ ∼ ⇓ ⇓ ⇓ ⇓		⇓ ⇓ ⇓ ⇓	∼ ⇓ ⇓ ∼ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓	↓ ↓ ↑ ↑	↓ ∼ ↑ ∼ ↑ ↑ ↓ ↓ ↓ ∼ ∼ ∼ ↓	↓ ↑ ↓ ↑	∼ ↑ ↓ ↓ ∼ ↑ ∼ ∼ ↓ ↑ ↓ ↓ ↑
	33. ⇓	∼ ∼ o r3 o r7 o r14 x CodY x TnrA ⇓ ↑ ∼ o r3 o r7 o r14 x CodY x TnrA
	13.2. Unsafe Changes for ∆ 6 34. ∼ 68. ⇓ 99. ⇓	⇓		⇓ ⇓	⇓			⇓ ⇓	∼		∼ ↑	↑	↑ ↓	∼
	35. ∼ 69. ⇓ 100. ⇓	⇓		⇓ ⇓	∼		⇓ ⇓	⇓		↑ ↑	↓	∼ ↑	∼
	x CodY x TnrA ↓ ∼ ∼ ⇓ ∼ o r1 o r3 o r7 o r14 x CodY x TnrA 1. 36. ∼ ⇓ ↓
	70. ⇓	2. o r3 o r7 o r14 x CodY x TnrA ∼ ↓ ⇓ ∼ ∼ ↓ ↑
	37. ⇓ 71. ⇓	3. ⇓	∼	∼ ⇓ ∼	∼	↑ ∼	↑	↑	↑
	38. ⇓ 72. ⇓	4. ∼	∼	⇓	↑	⇓	∼	∼ ↑	↓	∼	↓
	39. ∼ 73. ⇓	∼		⇓	⇓			⇓	∼		↓	↓	∼	↑
	40. ∼ 74. ⇓	∼	o r1 o r3 o r7 o r14 ⇓ ⇓ ∼ ⇓ ∼ ↑	↓	↓
	5. ⇓ 41. ∼ ⇓ 75. ⇓ ∼ ∼	∼ ∼ ⇓ ⇓	∼	⇓ ↓	↑	↓
	6. ∼ 42. ∼ ⇓ 76. ⇓ ∼ ∼	⇓ ⇓	⇓	∼	↑	⇓ ↓	∼	↑
	7. ∼ ∼ ∼ ∼ o r1 o r3 o r7 x CodY x TnrA ⇓ ⇓ ⇓ ↑ ↓ ∼ ∼ ⇓ ↑ ↑ o r1 o r3 o r7 o r14 x CodY x TnrA 77. ⇓ 78. ⇓ 8. ⇓ ∼ ∼ ∼ ↓ 43. ⇓ ∼ ∼ ↓ ↓ 79. ∼ ⇓ ⇓ ∼ ↑ ↓ ∼ 9. ⇓ ∼ ∼ ∼ ∼ 44. ⇓ ∼ ∼ ↓ ↑ ↓ 10. ⇓ ∼ ∼ ∼ ∼ 45. ⇓ ∼ ∼ ↑ ↓ o r3 o r7 o r14 x CodY x TnrA ↑ 11. ∼ ⇓ ∼ ∼ ∼ 46. ⇓ ∼ ∼ ↑ ↑ 80. ⇓ ⇓ ∼ ↑ ↑ ↑ 12. ∼ ⇓ ∼ ∼ ↑ 47. ∼ ⇓ ∼ ↓ ↑ 81. ⇓ ∼ ⇓ ↓ ↓ ∼ 13. ∼ ∼ ⇓ ∼ ↓ 48. ∼ ⇓ ∼ ↑ ↓ 82. ⇓ ∼ ⇓ ↓ ↑ ∼ 14. ∼ ∼ ⇓ ∼ ∼ 49. ∼ ⇓ ∼ ↑ ↑ 83. ⇓ ∼ ⇓ ↑ ↓ ↓ 15. ∼ ∼ ⇓ ∼ ↑ 50. ∼ ∼ ⇓ ↓ ↓ 84. ⇓ ∼ ⇓ ↑ ↑ ∼ 16. ∼ ∼ ∼ ⇓ ↓ 51. ∼ ∼ ⇓ ↓ ↑ 85. ∼ ⇓ ⇓ ↓ ↓ ∼ 17. ∼ ∼ ∼ ⇓ ∼ 52. ∼ ∼ ⇓ ↑ ↓ 86. ∼ ⇓ ⇓ ↓ ↑ ↓ 18. ⇓ ∼ 19. ⇓ ↑ 55. ∼ ⇓ ↓ ↑ ∼ 54. ∼ ⇓ ↓ ↓ ↑ 53. ⇓ ∼ ↑ ↑ o r14 x CodY x TnrA 87. ∼ ⇓ ⇓ ↑ ↓ o r7 o r14 x CodY x TnrA 88. ∼ ⇓ ⇓ ↑ ↑
	56. ∼ 20. 57. ∼	⇓ x CodY x TnrA ↑ ↓ ↓ ⇓ ↑	↓ ↑
			21.			↓			↑
			22.			↑			↓
			23.			↑			↑

  var Directions: x_TnrA; constraint x_TnrA in {Do, Up, No}; var Directions: o_1; constraint o_1 in {DoN, No, NoN}; var Directions: o_14; constraint o_14 in {DoN, No, NoN}; var Directions: o_3; constraint o_3 in {DoN, No, NoN}; var Directions: o_7; constraint o_7 in {DoN, No, NoN}; var Directions: o_8; constraint o_8 in {DoN, No, NoN};

	constraint let {
	%% local variables %%
	var Directions: y_Leu; constraint y_Leu in {Do, Up, No};
	var Directions: z_BSCodY;
	var Directions: z_PIlv_Leu;
	var Directions: z_BSCcpA;
	var Directions: z_BSTnrA;
	var Directions: z_OpBkLBcd;
	var Directions: w_d4e6;
	var Directions: w_d1e43;
	var Directions: w_d1e57;
	var Directions: w_d1e59;
	var Directions: w_d1e64;
	var Directions: w_d4e19;
	var Directions: w_d1e97;
	var Directions: w_d4e33;
	var Directions: w_d1e115;
	var Directions: w_d4e51;
	var Directions: w_d1e205;
	var Directions: w_d1e217;

Now BS CcpA may inhibit the acceleration of r

through CcpA via the new reaction r 9 , that replaces the previous direct acceleration of r 2 by CcpA.

Acknowledgements

The authors would like to thank Emilie Allart and Debarun Dhali for proofreading. This research was supported by the European Union (Marie Curie ITN AMBER, 317338), the French National Research Agency (ICEBERG-ANR-10-BINF-06-01), and the University of Lille, through the BQR project "Biologie synthétique pour la synthèse dirigée de peptides microbiens bioactifs".

15. MiniZinc Version of Difference Constraint over ∆ 6 int: No = 1; int: Up = 2; int: Do = 3; int: NoN = 4; int: UpN = 5; int: DoN = 6; set of int: Directions = {No, Up, Do, NoN, UpN, DoN}; array [1..6] of string: dirStrings = ["No", "Up", "Do", "NoN", "

UpN", "DoN"]; include "table.mzn"; % s1 (presumably a reaction flow) is accelerated by s2 ( presumably a chemical species) predicate abstractAccelerator(var int: s1, var int: s2, var int: finalRes) = abstractAcceleratorTable(s1, s2, finalRes); % abstractAcceleratorBackward(s1, s2, finalRes) /\ abstractAcceleratorForward(s1, s2, finalRes);