François Fages

Steven Gay

Sylvain Soliman

Inferring Reaction Systems from Ordinary Differential Equations

Keywords: systems biology, chemical reaction network theory, ordinary differential equations, SBML

In Mathematical Biology, many dynamical models of biochemical reaction systems are presented with Ordinary Differential Equations (ODE). Once kinetic parameter values are fixed, this simple mathematical formalism completely defines the dynamical behavior of a system of biochemical reactions and provides powerful tools for deterministic simulations, parameter sensitivity analysis, bifurcation analysis, etc. However, without requiring any information on the reaction kinetics and parameter values, various qualitative analyses can be performed using the structure of the reactions, provided the reactants, products and modifiers of each reaction are precisely defined. In order to apply these structural methods to parametric ODE models, we study a mathematical condition for expressing the consistency between the structure and the kinetics of a reaction, without restricting to Mass Action law kinetics. This condition, satisfied in particular by standard kinetic laws, entails a remarkable property of independence of the influence graph from the kinetics of the reactions. We derive from this study a heuristic algorithm which, given a system of ODEs as input, computes a system of reactions with the same ODE semantics, by inferring well-formed reactions whenever possible. We show how this strategy is capable of automatically curating the writing of ODE models in SBML, and present some statistics obtained on the model repository biomodels.net. 1 .

Introduction

In Mathematical Biology, many models are presented as a system of Ordinary Differential Equations (ODEs). Once the kinetic parameter values are fixed, this simple mathematical formalism completely defines the dynamical behavior of a system of biochemical reactions. It provides powerful tools for both transient and steady-state analysis via numerical integration, parameter sensitivity analysis, or bifurcation analysis, but only when kinetic information is available.

In absence of knowledge on the kinetics of each reaction, various qualitative analyses can nevertheless be performed using the structure of the reactions. This approach has rapidly developed in Systems Biology for reasoning on large interaction networks, with for instance, the analysis of qualitative attractors in a logical dynamics of gene networks à la Thomas [START_REF] Thomas | A complex control circuit : regulation of immunity in temperate bacteriophages[END_REF][START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata : a logical description[END_REF], reachability and temporal logic properties in reaction networks [START_REF] Eker | Pathway logic: Symbolic analysis of biological signaling[END_REF][START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF][START_REF] Chabrier-Rivier | Modeling and querying biochemical interaction networks[END_REF][START_REF] Fages | Modelling and querying interaction networks in the biochemical abstract machine BIOCHAM[END_REF][START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF], structural invariants in the Petri net representation of the reactions [START_REF] Reddy | Petri net representations in metabolic pathways[END_REF][START_REF] Zevedei-Oancea | Topological analysis of metabolic networks based on petri net theory[END_REF][START_REF] Angeli | A petri net approach to persistence analysis in chemical reaction networks[END_REF][START_REF] Chaouiya | Petri net modelling of biological regulatory networks[END_REF][START_REF] Rohr | Snoopy -a unifying petri net framework to investigate biomolecular networks[END_REF][START_REF] Soliman | Invariants and other structural properties of biochemical models as a constraint satisfaction problem[END_REF], or model reductions using graph theory concepts [START_REF] Gay | A graphical method for reducing and relating models in systems biology[END_REF][START_REF] Gay | On the subgraph epimorphism problem[END_REF]. These qualitative analysis tools do not rely on kinetic information, but on the structure of the reaction network which has thus to be correctly written as a set of formal reactions, with well-identified reactants, products and modifiers (and in certain cases their stoichiometry) for each reaction.

For instance, in [START_REF] Kaleta | Using chemical organization theory for model checking[END_REF], it is elaborated that structural information hidden in kinetic laws may affect the results obtained from structural analyses, such as elementary mode analysis [START_REF] Schuster | A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks[END_REF], flux balance analysis [START_REF] Varma | Metabolic flux balancing: basic concepts, scientific and practical use[END_REF], chemical organization theory [START_REF] Dittrich | Chemical organisation theory[END_REF], deficiency analysis or chemical reaction network theory [START_REF] Feinberg | Mathematical aspects of mass action kinetics[END_REF][START_REF] Shinar | Structural sources of robustness in biochemical reaction networks[END_REF].

It is worth noticing that these structural analyses may also directly support dynamic analyses. For instance, [START_REF] Koh | A decompositional approach to parameter estimation in pathway modeling: a case study of the akt and mapk pathways and their crosstalk[END_REF] applies network decomposition for a modular parameter estimation approach, [START_REF] Angeli | A petri net approach to persistence analysis in chemical reaction networks[END_REF] introduces a structural persistence criterion, Petri net place invariants reveal conservation laws in [START_REF] Soliman | Finding minimal P/T-invariants as a CSP[END_REF], while transition invariants can be used to identify fragile nodes and the core of a network [START_REF] Grafahrend-Belau | Modularization of biochemical networks based on a classification of petri net by T-invariants[END_REF], or to determine steady state solutions [START_REF] Nabli | Steady-state solution of biochemical systems, beyond S-Systems via T-invariants[END_REF].

Furthermore, knowing the correct structure of each reaction is mandatory when a reaction network must be interpreted as a stochastic process (Continuous-Time Markov Chain, CTMC) à la Gillespie [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF].

The question of the correct identification of a structured reaction model from a system of ODEs is thus important and is not new. Actually for the restricted case of models with only Mass Action kinetics a general solution is provided in [START_REF] Hárs | On the inverse problem of reaction kinetics[END_REF]. This approach was evolved over the years, see for instance [START_REF] Szederkényi | Inference of complex biological networks: distinguishability issues and optimization-based solutions[END_REF] for sparse/dense/core solutions when numerical values are provided for the parameters, or [START_REF] Soliman | A unique transformation from ordinary differential equations to reaction networks[END_REF] for unicity conditions in the symbolic case, still in the restricted framework of mass action law kinetics. In [START_REF] Kaleta | Using chemical organization theory for model checking[END_REF], the authors present an algorithm that uncovers hidden structural information for some Systems Biology Markup Language (SBML) [START_REF] Hucka | The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models[END_REF][START_REF] Hucka | Systems biology markup language (SBML) level 2: Structures and facilities for model definitions[END_REF] models of the biomodels.net repository [START_REF] Le Novère | BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF], with restricting to reaction models without inhibitors.

In this paper, we describe an algorithm for finding a reaction models for a given system of ODEs, considering reaction with inhibitors and general kinetic expressions. The first contribution of this paper is to propose a mathematical condition for expressing the consistency between the kinetic expression and the reactant-product-inhibitor structure of a reaction. We introduce well-formedness (Def. 2.5) and strictness (Def. 2.6) conditions for reactions, and show that they are satisfied by standard kinetics such as Mass Action law, Michaelis-Menten, Hill and negative Hill kinetics. The well-formedness condition is also shown to entail a property of independence of the influence graph (or symbolic Jacobian matrix) from the kinetics of the reactions (Thm. 2.16). This result generalizes a previous result in [START_REF] Fages | From reaction models to influence graphs and back: a theorem[END_REF][START_REF] Fages | Abstract interpretation and types for systems biology[END_REF] to reactions with inhibitors. It shows that the influence graph of a well-formed reaction system with inhibitors is essentially independent of the kinetics, can be computed in linear time in the number of reactions when the number of species per reaction is bounded, and can thus advantageously be used to perform multi-stationarity analyses by circuit analysis à la Thomas [START_REF] Soliman | A stronger necessary condition for the multistationarity of chemical reaction networks[END_REF][START_REF] Soulé | Mathematical approaches to differentiation and gene regulation[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata : a logical description[END_REF][START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF][START_REF] Thomas | A complex control circuit : regulation of immunity in temperate bacteriophages[END_REF].

The second contribution of this paper is to use these well-formedness and strictness conditions to prove the completeness of a new general algorithm for inferring a reaction system equivalent to an ODE system. This algorithm, of time complexity in O(n × t) where n is the number of variables and t the number of terms in the ODE, is shown to preserve the ODE semantics of the reactions (soundness 3.10), as well as their well-formedness when applied to the ODE semantics of a non-decomposable well-formed reaction system (weak completeness 3.12).

Our third contribution is to show that our algorithm can be used to automatically curate the writing of ODE models with reactions, as required in SBML. The fact that SBML has become a standard for sharing and publishing models has helped in making modelers formalize the reaction structure of their models. Unfortunately, SBML does not enforce any strong coherence between the structure and the kinetics of a reaction. Therefore the structural interpretation of models transcribed in reaction-based formalisms such as SBML may vary according to different choices of representation of the original ODE model as a reaction system, and may invalidate some structural analyses. We compare our results to the one presented in [START_REF] Kaleta | Using chemical organization theory for model checking[END_REF], and provide some statistics obtained on the rewriting in SBML of the curated part of the biomodels.net repository, showing that our method is able to automatically decrease the number of non well-formed reaction systems from 65% to 29%.

A Theory of Well-formed Reactions and Kinetics

In this section, we consider a finite set S = {x 1 , . . . , x s } of s molecular species, and a finite set of n reactions over S which are formally represented as multiset rewriting rules with kinetic expressions.

Multisets are used for representing reactants, products and inhibitors in reactions. A multiset s of molecular species is a function S -→ N which gives the number (stoichiometric coefficient) s(x) of each molecular species x ∈ S in s. We have s(x) = 0 if x does not belong to s, and s(x) ≥ 1 if x belongs to s, which is also written x ∈ s by abuse of notation. The empty multiset is written ∅. Equivalently, a multiset s will also be denoted by the linear expression m i=1 s(x i ) × x i , which gives the stoichiometric coefficients of each molecular species x i in s. This corresponds to the classical chemical notation

2H + O -→ H 2 O.
We shall now introduce the well-formedness and strictness conditions and describe some of their properties.

Well-formedness and Strictness Conditions

In the following definition, a reaction is composed of multisets for reactants, products and inhibitors that are not assumed to be disjoint. Definition 2.1. A reaction is a quadruple (r, m, p, f ), where r is the multiset of reactants, m the multiset of inhibitors, p the multiset of products, and f , called kinetic expression, a mathematical function over molecular species concentrations, f : R s -→ R. A reaction system is a finite set of reactions.

The species that are both reactants and products in a reaction are called catalysts. For the sake of readability, a reaction (r, m, p, f ) will also be written r / m f -→ p or just r f -→ p if it has no inhibitor, i.e. when m = ∅. The kinetic expression will also be omitted if it is not relevant.

Example 2.2. For instance, the following reaction, transcribed from Kohn's map of the cell cycle [START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF],

pMPF + Cdc25 k 1 ×pMPF ×Cdc25 -→ MPF + Cdc25
expresses the activation of the Mitosis Promoting Factor MPF by the kinase Cdc25 . It has as rate law f = k 1 × pMPF × Cdc25 , i.e. a Mass Action kinetics with parameter k 1 . In this reaction, pMPF is a reactant, MPF a product, Cdc25 a reactant and a product at the same time, i.e. a catalyst in our terminology, and there is no inhibitor.

A simplified version of that reaction can be written by omitting the kinase Cdc25 , as follows:

pMPF V ×pMPF /(K+pMPF )

-→ MPF

That form typically derives from three reactions describing the reversible association of pMPF and Cdc25 and the dissociation to MPF , by making a quasi steady state approximation on Cdc25 , which results in a Michaelis-Menten kinetics with parameters K and V.

It is worth noting that in a reaction, a reactant or a product can also be an inhibitor if it appears in m.

Example 2.3. For instance, the Botts-Morales general modifier mechanism accounts for a modifier M that can enhance and slow down a reaction A -→ B, depending on its concentration [START_REF] Katsumata | Graphic representation of botts-morales equation for enzyme-substrate-modifier system[END_REF]. This can be represented in our setting by a reaction of the form

A + M / M -→ B + M.
SBML does not distinguish between catalysts and inhibitors which are just considered as "modifiers" in SBML annotations. However we find it useful for the theory developed here to distinguish between the activation or inhibitory effects of a modifier, and mark it syntactically as such in the structure of the reaction. If a modifier has both activation and inhibitory effects, it will just appear in r, m and p in our setting, without loss of generality.

It is also worth noting that we consider only irreversible reactions, as in Feinberg's Chemical Reaction Network theory [START_REF] Feinberg | Mathematical aspects of mass action kinetics[END_REF]. A reversible reaction is thus represented by two reactions, one for each direction. This is one important difference with the Systems Biology Markup Language (SBML) that permits the declaration of a reversible reaction with only one single kinetic expression which can be negative.

These distinctions do not affect the system of ODEs that is classically associated to a reaction system by the Reaction Rate Equation as follows:

Definition 2.4. The ODE semantics of a reaction system

R = { r i / m i f i -→ p i } i=1,...,n
over molecules {x 1 , . . . , x s }, is the system of ordinary differential equations

ẋj = n i=1 (p i (x j ) -r i (x j )) × f i for 1 ≤ j ≤ s.
Our aim is to go in the reverse direction, that is to infer from any ODE system a reaction system with the same ODE semantics. Let us first remark that any ODE system ẋj = g j can be trivially transcribed in a reaction system using artificial synthesis reactions for each molecular species, with the terms of the differential equation as kinetic expressions, as follows:

∅ g j -→ x j
Since the ODE semantics is identical to the original ODE system, this is correct as far as numerical simulations are concerned, but prevents the use of structural analysis methods or stochastic simulations as the structures of the reactions are totally meaningless. It is worth remarking that some ODE models have nevertheless been transcribed in SBML using that scheme, since it does not affect simulations. This is the case for instance of model BIOMD0000000008.xml in biomodels.net for the ODE model of [START_REF] Gardner | A theory for controlling cell cycle dynamics using a reversibly binding inhibitor[END_REF]. We will use that example in Section 3.1 to illustrate our reaction inference algorithm and its capability of automatically curating the writing in SBML of ODE models as reaction systems.

In order to try to infer meaningful reactions from ODEs, we are interested in mathematical conditions for expressing the consistency of the kinetic expression f with the structure (r, m, p) of a reaction. Furthermore, since it is common practice to aggregate a system of elementary reactions in one abstract reaction with more complex kinetics (the simplest example of which are Michaelis-Menten and Hill kinetics for enzymatic reactions), we do not content ourselves with elementary kinetic expressions such as mass action law kinetics, but seek abstract consistency properties that can be applied to any mathematical expression given as kinetics. This is in contrast to most work on chemical reaction network theory [START_REF] Feinberg | Mathematical aspects of mass action kinetics[END_REF][START_REF] Shinar | Structural sources of robustness in biochemical reaction networks[END_REF][START_REF] Soliman | A unique transformation from ordinary differential equations to reaction networks[END_REF], but in accordance with the use in SBML of MathML for writing the kinetic expressions without any limitation on the use of mathematical symbols. Definition 2.5. A reaction (r, m, p, f ) over molecular species {x 1 , . . . , x s } is well-formed if the following conditions hold:

1. f (x 1 , . . . , x s ) is a partially differentiable function, non-negative on R s + ; 2. x i ∈ r if and only if ∂f /∂x i ( x) > 0 for some value x ∈ R s + ; 3. x i ∈ m if and only if ∂f /∂x i ( x) < 0 for some value x ∈ R s + . The first condition expresses that the kinetic expression must be a differentiable and non-negative function for all non-negative values of the variables. The second (resp. third) condition states that the partial derivative of f w.r.t. a reactant (resp. an inhibitor) must be positive (resp. negative) for some (not necessarily all) non-negative values of the variables.

It is worth noting that we do not impose the monotonicity condition that for any variable x i ∈ V, ∂f /∂x i should be either non-negative on the positive orthant, or non-positive on the positive orthant. In our setting, a molecular species can thus be both a reactant and an inhibitor in a well-formed reaction, depending on the values of the concentrations. On the other hand we shall make use of the following: Definition 2.6. A reaction (r, m, p, f ) is strict if its kinetics f (x 1 , . . . , x s ) = 0 whenever x j = 0 for any x j such that r(x j ) > 0.

This condition expresses that the kinetics must be zero if the concentration of one of the reactants is zero. If the kinetics is a rational expression, that strictness condition implies that the kinetic expression is a product of the reactants with a fractional expression defined for all non-negative values of the variables. More generally it enforces the positivity of the system:

Definition 2.7 (Positive System). A dynamical system over R k is called positive if R k
+ is an invariant set for the system, i.e., ∀x 0 ≥ 0, t ≥ 0 x(t, x 0 ) ≥ 0 Proposition 2.8 (Positivity). The ODE semantics of a well-formed and strict reaction system defines a positive system.

Proof. In Definition 2.4 we have ẋj = n i=1 (p i (x j ) -r i (x j )) × f i and since the system is well-formed, the f i are all non-negative. The only negative terms thus have r i (x j ) > 0 and from the strictness condition this entails that f i = 0 when x j = 0. Hence ẋj ≥ 0 whenever x j = 0 since it is a sum of non-negative terms. Therefore x j cannot become negative when its initial value is non-negative, and since this holds for all j, the system is positive.

The strictness condition excludes the writing of a reversible reaction with one single reaction by summing the kinetic expressions of each direction, as allowed in SBML, when the reactants differ from the products.

It also excludes the existence of a strict well-formed reaction system for any ODE system, as shown by Example 2.9. The equation ẋ = -k is not the ODE semantics of any strict well-formed reaction system, since that ODE defines a non-positive system (Prop. 2.8). That ODE can be associated to the non-strict well-formed reaction system

x l×x -→ 2 × x x k+l×x -→ ∅
(where the kinetic expression is not null when x = 0). This is the result computed in that case by Algorithm 3.6 described later.

Example 2.10. As for an example with inhibitors, let us consider the following three reactions representing the core of the action of the Circadian clock on the Cell Cycle, as described by Matsuo et al. [START_REF] Matsuo | Control mechanism of the circadian clock for timing of cell division in vivo[END_REF]:

pMPF + Cdc25 k 1 ×pMPF ×Cdc25 -→ MPF + Cdc25 MPF + Wee1 k 2 ×MPF ×Wee1 -→ pMPF + Wee1 ∅/Clock k 3 /(k 4 +Clock ) -→ Wee1
where k 1 , k 2 , k 3 are parameters. The first reaction is the one of Example 2.2. Those reactions are well-formed and strict. In particular, we have

∂ Ẇee1 /∂Clock = ∂ k 3 /(k 4 + Clock ) /∂Clock = -k 3 /(k 4 + Clock ) 2 < 0 for showing the inhibitory effect of Clock in the synthesis reaction of Wee1. Their ODE semantics is ṗMPF = k 2 × MPF × Wee1 -k 1 × pMPF × Cdc25 ṀPF = k 1 × pMPF × Cdc25 -k 2 × MPF × Wee1 Ẇee1 = k 3 /(k 4 + Clock ) Ċdc25 = 0 Ċlock = 0
The well-formedness and strictness conditions are satisfied by standard kinetic laws. One can easily check Proposition 2.11. Reactions with mass action law kinetics:

j n j × x j k× j x n j j -→ p Michaelis-Menten kinetics: x V ×x/(K+x) -→ y Hill kinetics: x V ×x n /(K n +x n ) -→ y
or negative Hill kinetics:

∅/x V /K n +x n -→ y
with rate constants k, V, K > 0 and exponent n ≥ 1, are well-formed and strict.

We shall see in Section 4.2 that these conditions are currently violated in a majority of reaction systems of the biomodels.net repository, but that most of them can be automatically corrected by modifying their structure and writing in SBML, without changing their ODE semantics.

Influence Graph associated to a Well-formed Reaction System

The influence graph between molecular species induced by the ODE semantics of a wellformed reaction system enjoys a remarkable property of independence from the kinetics, which we present in this section. Influence graphs have been initially introduced in the setting of gene regulatory networks [START_REF] Thomas | A complex control circuit : regulation of immunity in temperate bacteriophages[END_REF] as a simple abstraction enabling reasoning about complex regulation mechanisms. These graphs completely abstract from the precise interactions, especially at post-transcriptional level, and retain only the activation and inhibitory effects on gene transcription. As conjectured in [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF], the existence of a positive circuit (resp. a negative circuit) in an influence graph has been proved to be a necessary condition for multi-stationarity, e.g. for cell differentiation, (resp. for oscillations, e.g. for homeostasis) in different formalisms, and in particular for ODE systems in [START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF][START_REF] Soulé | Mathematical approaches to differentiation and gene regulation[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF] and recently in [START_REF] Soliman | A stronger necessary condition for the multistationarity of chemical reaction networks[END_REF] for the ODE semantics of non-linear reaction systems.

Here, we show that in a well-formed reaction system, and under a very general assumption, the influence graph of the reactions is identical to the influence graph of the ODE semantics of the reactions.

On the one hand, in an ODE system, the influence graph is mathematically defined by the signs of the coefficients in the Jacobian matrix of the system, (∂ ẋi /∂x j ), as follows: Definition 2.12. The differential influence graph (DIG) associated to (the ODE semantics of ) a reaction system is the graph that has for vertices the molecular species, and for labeled edges the following set of signed edges: where the positive influences are represented by green arrows with triangular tips, and negative influences are represented by red arrows with blunt tips. For instance, the negative influence of Clock on Wee1 comes from the negative sign of ∂ Ẇee1 /∂Clock as detailed in 2.10. There are negative loops on MPF and MPF since ∂ ṀPF /∂MPF = -k 2 Wee1 < 0 and ∂ ṗMPF /∂pMPF = -k 2 Cdc25 < 0, and not on Cdc25 since Ċdc25 = 0. Note that a useful circuit analysis in this example would necessitate considering the reactions of formation of MPF and is beyond the scope of this paper.

{x i -→ + x j | ∂ ẋj /∂x i ( x) > 0 for some value x ∈ R s + } ∪ {x -→ -x j | ∂ ẋj /∂x i ( x) < 0 for some value x ∈ R s + } Example 2.
On the other hand, in a reaction system, one can define an influence graph directly from the stoichiometry of the reactions, ignoring the kinetics, as follows:

Definition 2.14. The stoichiometric influence graph (SIG) associated to a finite set R of reactions is the graph that has for vertices the molecular species, and for labeled edges the following set of signed edges:

{x -→ + y | either p i (y) -r i (y) > 0 and r i (x) > 0, or p i (y) -r i (y) < 0 and m i (x) > 0, for some reaction i} ∪ {x -→ -y | either p i (y) -r i (y) < 0 and r i (x) > 0, or p i (y) -r i (y) > 0 and m i (x) > 0, for some reaction i}

Intuitively, there is a positive (resp. negative) arc from x to y if x is a reactant in a reaction that produces more (resp. less) y than it consumes, or an inhibitor in a reaction that consumes more (resp. less) y than it produces.

Unlike the DIG, which needs to compute the sign of partial derivatives, the SIG can be easily computed in linear time in the number of reactions, assuming that the number of species per reaction is bounded, since it is sufficient to parse the stoichiometric coefficients of the reactions. As already shown in [START_REF] Fages | From reaction models to influence graphs and back: a theorem[END_REF], the SIG is an over-approximation of the DIG: Theorem 2.15 ([36]). For any finite set R of well-formed reactions, the DIG of R is a subgraph of the SIG of R.

We show here that, even in the presence of inhibitors, the SIG is in fact identical to the DIG with an extra assumption. Let us say that a tuple of molecular species (x, y) is in conflict in an influence graph if we have both x -→ + y and x -→ -y. Theorem 2.16. For any finite set R of well-formed reactions such that the SIG of R contains no conflict, the DIG and the SIG are identical.

Proof. We just have to prove that the SIG is a subgraph of the DIG. Let us consider an arc x -→ + y in the SIG. By Definition 2.14 there exists a reaction i with either p i (y) -r i (y) > 0 and r i (x) > 0, or p i (y) -r i (y) < 0 and m i (x) > 0. Since the reaction is well-formed, we have either p i (y) -r i (y) > 0 and ∂f i /∂x( z) > 0, or p i (y) -r i (y) < 0 and ∂f i /∂x( z) < 0, for some z ∈ R s + . Now, if p i (y)-r i (y) > 0 then f i occurs in ẏ with a positive sign. Since ∂f i /∂x( z) > 0 and there is no conflict in the SIG, we thus get ∂ ẏ/∂x( z) > 0, i.e. x -→ + y is in the DIG. Similarly, if p i (y) -r i (y) < 0, f i occurs in ẏ with a negative sign and ∂f i /∂x( z) < 0, hence ∂ ẏ/∂x( z) > 0, i.e. x -→ + y is in the DIG. The proof for an arc x -→ -y in the SIG is symmetrical.

Corollary 2.17. The DIG of a finite set of well-formed reactions without conflict in its SIG, is independent of the kinetic expressions.

Corollary 2.18. The DIG of a finite set of well-formed reactions without conflict in its SIG, is computable in linear time in the number of reactions, when the number of species appearing in a reaction is bounded.

The SIG of Example 2.10 is trivial to compute and since it contains no conflict, we can predict by Theorem 2.16 that it is identical to its DIG depicted in Example 2.13.

Example 2.19. As for an example of conflict, in the simplified model of the yeast cell cycle of [START_REF] Tyson | Modeling the cell division cycle: cdc2 and cyclin interactions[END_REF], the double activation reactions of MPF through Cdc25 and Wee1 [START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF], are simplified in a single autocatalytic reaction in parallel with a deactivation reaction:

pM P F + M P F -→ 2 × M P F M P F -→ pM P F
Such reactions create a conflict in the SIG, namely MPF -→ -pMPF and MPF -→ + pMPF . In general, there is a possibility that such conflicting direct influences in the SIG may be balanced in the ODEs and do not appear in the DIG. This situation is however quite pathological and rare in practice, and occurs when over-simplifications are made. For instance, Kohn's map of the cell cycle control [START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF] contains 800 reactions [START_REF] Chabrier-Rivier | Modeling and querying biochemical interaction networks[END_REF] and does not contain any conflict in its SIG [START_REF] Fages | From reaction models to influence graphs and back: a theorem[END_REF]. The conflict of influences between MPF and pMPF in Tyson's model comes from the compression in one loop of the two positive circuits through Wee1 and Cdc25 respectively. The decompression of this loop makes disappear the influence conflict.

Thomas's necessary condition for a system to exhibit multi-stationarity is the existence of a positive circuit, i.e., a simple oriented cycle such that the product of the signs of its edges is positive, in the DIG [START_REF] Soulé | Graphic requirements for multistationarity[END_REF]. That condition has proven useful to reason about gene interaction networks and predict the possibilities of multi-stationarity, i.e. cell differentiation. However, Thomas's original condition provides no information in presence of reactions with two reactants, since a reaction like for instance A + B -→ C immediately creates a positive circuit of negative influences between A and B in the associated SIG and DIG for any reasonable kinetics. This counter-example has been recently rule out in [START_REF] Soliman | A stronger necessary condition for the multistationarity of chemical reaction networks[END_REF], where it is shown that Thomas's conditions can be made stronger for reactions models, by labeling the influence edges by the reactions they come from, and by restricting the analysis of circuits to circuits labeled by different reactions. With this stronger condition for multi-stationarity, the analysis of labeled circuits in the DIG of a reaction system does provide information on its capabilities of exhibiting multi-stationarity. Theorem 2.16 shows, perhaps surprisingly, that for well-formed reaction systems without conflicts, the DIG is essentially independent from the kinetics, and in fact identical to the SIG, which is easy to compute and can be used to perform multi-stationarity analysis by circuit analysis.

Reaction System Inference Algorithm

In this section we present an algorithm to infer a reaction system from an arbitrary ODE system, and study its properties. The algorithm proceeds in two steps: one first step for infering hidden molecules corresponding to linear invariants of the ODE system, and one second step for infering the reactions.

Motivating Example

As remarked in Section 2.1, any ODE model can be transcribed in a reaction system using artificial synthesis and degradation reactions for each molecular species, with the positive, respectively negative, terms of the differential equation for the variables as kinetic expressions. While preserving the ODE semantics and thus ODE simulations, such a transcription prevents the use of structural methods and stochastic simulations to analyze the system.

Such a transcription has nevertheless been used in biomodels.net to write the ODE model of [START_REF] Gardner | A theory for controlling cell cycle dynamics using a reversibly binding inhibitor[END_REF] in SBML and create BIOMD0000000008.xml. This model adds a control mechanism to the cell-cycle model of Goldbeter et al. in [START_REF] Goldbeter | A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase[END_REF] but with this transcription in SBML, the reaction graph is not even connected.

Here are some of the reactions of this model (after expansion of the macros used in the original writing) which illustrate the problem:

∅ (1-M )×C×V 1 ×(C+K -1 6 )/(K 1 +1-M ) -→ M M M ×V 2 /(K 2 +M ) -→ ∅ ∅ M ×V 3 ×(1-X)/(K 3 +1-X) -→ X X V 4 ×X/(K 4 +X) -→ ∅
One can notice that ∂f 1 /∂C = 0, where f 1 is the kinetic expression of the first reaction, but C is not a reactant nor an inhibitor. The model is therefore not well-formed.

One can also note that, though encoded in complicated MathML expressions, 1 -M (resp. 1 -X) appears in the synthesis of M (resp. X) as a way to represent the inactive form of M (resp. X). Indeed, [START_REF] Goldbeter | A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase[END_REF] states that "(1 -M) thus represents the fraction of inactive (i.e., phosphorylated) cdc2 kinase, while (1 -X) represents the fraction of inactive (i.e., dephosphorylated) cyclin protease".

When applied to the ODE system associated to this model, the reaction system inference algorithm presented in the next two sections, infers two hidden molecules and the following well-formed and strict reactions:

M i + C C×V 1 ×(C+1/K 6 )×M i /(K 1 +M i ) -→ M + C M V 2 ×M/(K 2 +M ) -→ M i X i + M M ×V 3 ×X i /(K 3 +X i ) -→ X + M X V 4 ×X/(K 4 +X) -→ X i
The two inactive forms are now explicitly represented by two inferred molecules, written M i and X i , and the actions of C on M and of M on X are properly transcribed. The reaction system inferred automatically from the ODE semantics is well-formed and strict, and in fact consistent with the graphical representation of the paper [START_REF] Goldbeter | A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase[END_REF] where dashed arrows represent catalytic effects: Proc. Natl. Acad. Sci Cyclin is synthesized at a constant rate (vi) and triggers the transformation of inactive (MI) into active (M) cdc2 kinase by enhancing the rate of a phosphatase (E1); a kinase (E2) reverts this modification. In the second cycle of the phosphorylation-dephosphorylation cascade, cdc2 kinase (identical to E3) elicits the transition from the inactive (X+) into the active (X) form of a protease that degrades cyclin; the activation of cyclin protease is reverted by a phosphatase (E4). Vi (i = 1-4) denotes the effective maximum rate of each of the four converter enzymes; vd denotes the maximum rate of cyclin degradation by protease X. As shown in Fig. 3, this minimal cascade is capable of autonomous oscillatory behavior.

keep the model simple and to allow for the straightforward generation of thresholds (see below), the formation of a complex between cyclin and cdc2 kinase will not be taken into account; instead, it is assumed that cyclin drives cdc2 activation by enhancing the velocity of an "activase" which (see the above discussion) might primarily represent a tyrosine (and, possibly, threonine) phosphatase. Such a direct activation of the phosphatase acting on phosphorylated cdc2 kinase is one of the hypothetical mechanisms originally put forward for cyclin action [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF][START_REF] Dittrich | Chemical organisation theory[END_REF]. A further assumption is that the maximum activity ofthe kinase inactivating cdc2-the cdc2 "inactivase" (7)-remains constant throughout the cell cycle.

That okadaic acid, an inhibitor of phosphatase 2A, behaves as a mitotic inducer has suggested that the phosphatase acting on cdc2 might be activated through phosphorylation and inac- tivated by phosphatase 2A (23-26). This minimal model will not take into account the possible modification of the activase, nor will it differentiate the roles of cyclins A and B, which appear ification is governed by the foll equations:

dC C = VI -VdX - dt i Kd + C dM (1-M) dt K1 + (1 -M) dX (1 -X) V =3 dt -K3 +(1-X) - C c
In the above equations, C denot tion, while M and X represent th kinase and the fraction ofactive cyc represents the fraction of inactive (i kinase, while (1 -X) represents th dephosphorylated) cyclin protease. Vd denote, respectively, the constan and the maximum rate of cyclin de reached for X = 1; Kd and & denot for cyclin degradation and for cycl phatase acting on the phosphorylat represents an apparent first-order nonspecific degradation of cyclin whose contribution is much smal degradation by protease X, is not n sole effect is to prevent the boundl conditions where the specific prot

The normalized parameters Vi an the kinetics of the enzymes E, (i = cycles of post-translational modific phosphatase (E1) and the kinase molecule, and on the other hand, th phosphatase (E4) acting on the cyc For each converter enzyme, the two the effective maximum rate and divided by the total amount of rel MT (total amount of cdc2 kinase) fo 9108 Cell Biology: Goldbeter

In that form, the inferred model is thus suitable for further structural analysis.

The following sections present the reaction system inference algorithm in two steps: first the algorithm for inferring hidden molecules corresponding, as above, to invariants, second the algorithm for inferring well-formed reactions whenever possible.

Inference Algorithm for Hidden Molecules

ODE models often contain algebraic invariants, i.e., algebraic equations relating variables of the model and that hold true in any solution of the ODE system. Among those, linear invariants λ i x i = Λ, e.g. mass conservation invariants, or Petri-net place invariants, are an important particular case. A linear invariant can be used to simplify a model by eliminating one variable and replacing it with a linear expression. This may have several advantages, but when writing the model with reactions, such simplifications performed on the ODE system need be reversed in order to restore the correct structure of the reactions on eliminated molecular species, as shown for instance in the previous section with the inactive forms M i of M , and X i of X.

A preprocessor is first applied before the reaction inference algorithm, in order to reverse the elimination of linear invariants and infer hidden molecules. The expressions f for which new molecules are introduced need be chosen with care in order to avoid the introduction of useless variables. Restricting the search to expressions of the form k -x or k -x -y where k is a constant or parameter, and x and y are molecule concentrations, has proven useful in practice. This leads to together with functional dependency equations z j = f j (x 1 , . . . , x s ) for 1 ≤ j ≤ k.

Proposition 3.2 (Soundness). Let O be an ODE system over variables {x 1 , . . . , x s }. The ODEs computed by Algorithm 3.1 for the time derivatives of x 1 , . . . , x s , are mathematically equivalent to the equations in O given the functional dependency equations z j = f j (x 1 , . . . , x s ) for the hidden molecules.

Proof. We prove that each step of the algorithm replaces equal by equal, and thus that the whole execution preserves the mathematical equivalence of the equations. First, step (1) is a purely syntactical transformation that does not change the ODE system O. Now note that all the other changes are of two forms. Either the introduction of a new variable z such that ż = λ i ẋi , together with the functional dependency equation z = k + λ i x i , steps (2a) and (3a). Since the differential equation on z is indeed the time derivative of the definition of z, this does not change the equations on ẋi . Or the replacement of k + λ i x i by z, steps (2b-d) and (3b-c), which are equal from the previous definition of z.

Inference Algorithm for Reactions

The inference algorithm for reactions is based on a syntactical normal form for ODE systems which facilitates the recognition of common subterms in the equations.

We consider ODEs and kinetic laws written in MathML as terms with mathematical operations and functions (e.g.+, -, /, ×, etc.), constants of R and variables representing species concentrations and parameters. It is beyond the scope of this paper to precisely describe the mathematical expressions allowed and the symbolic computation performed. However, let us call non-decomposable a term that:

• its functor (top function symbol) is neither + nor -;

• cannot be reduced at top-level by the algebraic laws of distributivity of the product and division on addition and subtraction, e.g. if its functor is × (resp. /) then the arguments (resp. the numerator) are not sums.

Definition 3.3. A reaction r / m f -→ p over molecular species {x 1 , . . . , x s } is non-decomposable if f is syntactically a non-decomposable term. Definition 3.4. A mathematical expression is in additive normal form if it is of the form k i=1 c i × t i
where c i are integers and t i are distinct non-decomposable terms without integer coefficients.

An ODE system is in additive normal form if each equation is in additive normal form, i.e. if it is of the form

ẋi = l j=1 c i,j × t j , 1 ≤ i ≤ s
where l is the number of non-decomposable terms t j in the system.

Additive normal forms are not unique, but any ODE system can be written in additive normal form through standard algebraic transformations (such as the distributivity of × over +). The non-decomposability condition excludes the composition of several reactions in a single one with a sum as kinetic expression. In particular, we have: Proposition 3.5. Any non-decomposable well-founded reaction system, such that its ODE semantics is a polynomial ODE system, is strict.

Proof. First notice that a polynomial kinetics once in additive normal form results in a sum of monomials as non-decomposable terms. Now, from the second condition of well-formedness in Definition 2.5, for each reaction (r, m, p, f ) we have r(x j ) > 0 implies ∃ x, ∂f /∂x j ( x) > 0, but since f is a monomial, this implies that f has degree at least 1 in x j , and therefore that f (x 1 , . . . , x s ) = 0 when x j = 0, i.e., (r, m, p, f ) is strict. Now, given an ODE system in additive normal form, the following algorithm can infer an equivalent reaction system by sorting the terms of the ODEs, and creating one reaction for each term (formalized in Prop. 3.9 below). This algorithms requires checking the sign of a partial derivative, and as described in Section 4.1, such checks can be arbitrarily difficult for arbitrary mathematical expressions, but can be over-approximated. We thus assume given a test program (exact or not) for testing the sign of partial derivatives: partial_has_pos_val(f, x) that answers if yes or no the partial derivative of the function f with respect to variable x takes a non-negative value for some input in R s , and such that ∃ y ∂f /∂x( y) > 0 ⇒ partial_has_pos_val(f, x). For computability reasons, the reverse implication is not required. These tests are used in Steps 4(c) and 4(d). ṗMPF + ṀPF = 0. One variable, e.g. pMPF , can thus be eliminated and replaced by c -MPF . This yields the following ODE system, where all k i are positive:

ṀPF = k 1 × (c -MPF ) × Cdc25 -k 2 × MPF × Wee1 Ẇee1 = k 3 /(k 4 + Clock ) Ċdc25 = 0 Ċlock = 0
When applied to this system, using the test for partial derivatives described in Section 4.1, Algorithm 3.6 infers the following reactions:

Cdc25 c×k 1 ×Cdc25 -→ Cdc25 + M P F M P F + Cdc25 k 1 ×Cdc25×M P F -→ Cdc25 M P F + W ee1 k 2 ×M P F ×W ee1 -→ W ee1 ∅/Clock k 3 /(k 4 +Clock) -→ W ee1
However, by applying first the hidden molecule inference Algorithm 3.1, a hidden molecular species M P F i is introduced for the expression c -M P F . This hidden molecule corresponds to the linear invariant M P F i + M P F = c. We have

ṀPF i = -k 1 × MPF i × Cdc25 + k 2 × MPF × Wee1
and when applied to this ODE system after the preprocessing step, Algorithm 3.6 now computes the correct reactions:

M P F i + Cdc25 k 1 ×M P F i×Cdc25 -→ M P F + Cdc25 M P F + W ee1 k 2 ×M P F ×W ee1 -→ M P F i + W ee1 ∅/Clock k 3 /(k 4 +Clock) -→ W ee1
By counting the loops, one can easily check Proposition 3.8 (Time complexity). On an ODE system O in additive normal form, Algorithm 3.6 computes a reaction system in time O(n × t), where n is the number of variables and t is the number of non-decomposable terms in O.

By executing symbolically the algorithm, one can similarly check that the result is characterized in mathematical terms by Proposition 3.9 (Inferred reactions). Given an ODE system in additive normal form with appearing terms T = {f 1 , . . . , f t }: ẋi = t u=1 c i,u × f u for 1 ≤ i ≤ s. The reaction system inferred by Algorithm 3.6 is the set of non-decomposable reactions

{ r u /m u fu -→ p u } 1≤u≤t where r u = {i | c i,u <0} (-c i,u ) × x i + {i | c i,u ≥0, partial_has_pos_val(fu,x i )} x i , p u = {i | c i,u >0} c i,u × x i + {i | c i,u ≥0, partial_has_pos_val(fu,x i )} x i
and m u is the set of variables x such that partial_has_pos_val(-f u , x). Theorem 3.10 (Soundness). The ODE semantics of the reaction system inferred by Algorithm 3.6 from an ODE system O is equal to O.

Proof. Let us suppose without loss of generality that O = { ẋi = t u=1 c i,u ×f u | 1 ≤ i ≤ s} is in additive normal form. The inferred reaction system is the set {r u /m u fu -→ p u } 1≤u≤t where

r u = {i | c i,u <0} (-c i,u ) × x i + {i | c i,u ≥0, partial_has_pos_val(fu,x i )} x i , p u = {i | c i,u >0} c i,u × x i + {i | c i,u ≥0, partial_has_pos_val(fu,x i )} x i ,
and m u is the set of variables y such that partial_has_pos_val(-f u , x i ).

The ODE system associated to these reactions is thus

{ ẋi = t u=1 (p u (x i ) -r u (x i )) × f u } 1≤i≤s = { ẋi = t u=1 c i,u × f u } 1≤i≤s = O.
Note that it does not depend on the test partial_has_pos_val. Algorithm 3.6 always computes a non-decomposable reaction system with an equivalent associated ODE system but this reaction system may not be well-formed. In particular, step 3b adds a variable x to the reactants of the reactions even if x does not appear in the kinetic expression f of the reaction. Therefore the algorithm may infer reactions with reactants that do not occur in the kinetic expression, as required for instance by Example 2.9.

We can measure the completeness of the method by showing that, at least, if we start from a well-formed reaction model, generate the ODE semantics, and from the ODE system solely, infer back a reaction model, the algorithm does infer a well-formed reaction model.

First, it is clear that the algorithm infers non-decomposable kinetics (Prop. 3.9) in which any variable appearing in the kinetics appears in the reaction as either reactant (step 4b, or step 4c for catalysts), inhibitor (step 4d) or both: Proposition 3.11. The reactions inferred by Algorithm 3.6 contain no reaction with a molecular species x appearing in the kinetic expression f with ∂f /∂x = 0, and not appearing as a reactant or inhibitor.

This proposition remains true even if the sets of variables for which the partial derivatives are positive (4c) or negative (4d) are over-approximated. However, for completeness an exact test is necessary. Theorem 3.12 (Weak completeness). When applied to the ODE semantics of a nondecomposable well-formed reaction system such that ∀1 ≤ i ≤ n, 1 ≤ j ≤ s ∂f i /∂x j > 0 ⇔ partial_has_pos_val(f i , x j ), Algorithm 3.6 does infer a non-decomposable well-formed reaction system. Furthermore, if the ODE system is polynomial, the inferred model is strict.

Proof. Let us consider the ODEs associated to a well-formed non-decomposable reaction

system R = { r i / m i f i -→ p i } i=1,...,n . The ODE system is of the form O = { ẋj = n i=1 (p i (x j ) -r i (x j )) × f i | 1 ≤ i ≤ m}
which is an additive normal form after evaluation of the integers p i (x j )-r i (x j ). By Prop. 3.9, the inferred reaction system is { r i /m i f i -→ p i } 1≤i≤n where f i is non-decomposable by hypothesis, r i = {j|p i (x j )<r i (x j )} (r i (x j ) -p i (x j )) × x j + {j | p i (x j )≥r i (x j ), ∂f i /∂x j >0} x j p i = {j|p i (x j )>r i (x j )} (p i (x j ) -r i (x j )) × x j + {j | p i (x j )≥r i (x j ), ∂f i /∂x j >0} x j , and m i = m i . Now for any variable x j , we have x j ∈ r i if and only if x j ∈ r i since either p i (x j ) < r i (x j ) or ∂f i /∂x j > 0. Similarly x j ∈ p i if and only if x j ∈ p i since either p i (x j ) > r i (x j ) or p i (x j ) = r i (x j ) and ∂f i /∂x j > 0, These equalities between the sets (not multisets) of reactants, products and inhibitors suffice to show the well-formedness of the inferred reactions.

Strictness in the polynomial case follows from Prop. 3.5.

Since we do not restrict ourselves to Mass Action kinetics, our algorithm may well infer reactions with other kinetic expressions in cases where purely Mass Action reactions were possible. This is an important difference between our algorithm and the previous algorithms, which are restricted to Mass Action kinetics [START_REF] Hárs | On the inverse problem of reaction kinetics[END_REF][START_REF] Szederkényi | Inference of complex biological networks: distinguishability issues and optimization-based solutions[END_REF][START_REF] Soliman | A unique transformation from ordinary differential equations to reaction networks[END_REF]. Furthermore, even if we restrict to polynomial ODEs and Mass Action kinetics for reaction, further conditions are necessary to grant the unicity of the solution [START_REF] Soliman | A unique transformation from ordinary differential equations to reaction networks[END_REF]. Furthermore, another Mass Action reaction system exists for this ODE system:

x k×x -→ x + 2 × y x k×x -→ ∅

Evaluation Results on biomodels.net

The ability to infer a reaction system from ODEs can be turned into some automatic curation algorithm, as was done in Theorem 3.12, by inferring the reactions from the ODE semantics of a starting reaction system. In this section, we evaluate this form of curation on repository of structured models.

Computability Issues

Since, like in SBML, we allow arbitrary mathematical expression for kinetic expressions, checking the well-formed conditions may raise arbitrary difficult symbolic computation problems. These conditions can be checked however by doing some approximations.

In our implementation in Biocham [1, 2]2 , the partial_has_pos_val proceeds as follows: the kinetic expressions are first normalized as if they were polynomials, stopping when a non-polynomial operator (anything else than +, -and ×) is found. For the polynomials, the exact computation of the sign of any partial derivative is easy. For the other terms, either they are recognized as a standard kinetics (like Hill functions) and once again the exact sign is extracted, or they are considered unknown and for any variable appearing we will assume that it is possible that ∂f /∂x becomes positive for some values, and negative for some values. This is a conservative over-approximation.

With these provisions, different syntactical conditions may indicate that a reaction is not well-formed. The conditions for a reaction to be ill-formed can be classified into three categories:

1. "K not R" indicates that the concentration of a compound appears in the kinetic law of a reaction, but this compound is neither a reactant nor an inhibitor of the reaction; 2. "R not K" indicates that some compound is marked as reactant or inhibitor in a reaction, but does not appear in the kinetic expression; 3. "Negative" indicates that a kinetic expression may be negative with non-negative concentration values.

Indeed, in a well-formed reaction with kinetic expression f , if a species x is neither a reactant nor an inhibitor, then ∂f /∂x = 0, hence x should not appear in the kinetic expression f . Similarly, if a species is a reactant or an inhibitor, then ∂f /∂x = 0, so x should appear in f . Moreover, f should be non-negative. These ill-formedness conditions are checked in Biocham using the previous approximations. They correspond to the warning messages that Biocham can raise when loading a reaction system.

Global analysis

The 424 models from the curated branch of the latest version (release 24) of the biomodels. net repository [START_REF] Le Novère | BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF] were used as benchmark to test our reaction system inference algorithm, and compare the results with the original writing of the models in SBML. Out of those 424 models only 361 define reactions with proper kineticLaws. The other ones only describe systems through events and rules, or with no kinetic information, and thus have no ODE semantics.

Our curation algorithm reads the SBML model, extracts the corresponding ODE system and infers from it a new reaction system. Table 4.2 summarizes the result of the procedure, as detected by Biocham warnings. Over the 361 reaction systems of the original curated part of biomodels.net with ODE semantics, our algorithm reveals hidden molecules in 58 models, 173 models with "K not R" warning, 123 models with "R not K" warning and 157 models with "negative kinetics" warning. Our algorithm is able to automatically curate the writing of these models with reactions by reducing the number of non well-formed models with a warning by more than the half, from 65% to 29%. Table 1: Number of models having a "K not R", "R not K", or "negative kinetics" warning among the original 361 models of the curated part of biomodels.net, and among the reaction systems automatically inferred from their ODE semantics. "Any warning" reflects model for which there was at least one of the three warning.

As predicted by Proposition 3.11, the Algorithm 3.6 completely removes the "K not R" warnings. For the two other warnings, since the algorithm focuses on non-decomposable kinetics, it results in curated models quite close to the original ones, but does not tackle thoroughly the case of reactions with rates independent of some reactant, for the reasons illustrated in Example 2.9 of for any other reason. Therefore, 103 over 361 models remain with a non well-formedness warning. [START_REF] Kaleta | Using chemical organization theory for model checking[END_REF] In [START_REF] Kaleta | Using chemical organization theory for model checking[END_REF], the authors also scan the biomodels.net repository and report finding 5 inconsistencies: models 44, 93, 94, 143 and 151. Their diagnostics is as follows, some reaction fluxes become negative during the simulations of those models because of missing reversibility indications in models 93, 94 and 143. In the two first cases they report that adding the reverse reactions makes the models consistent, whereas for 143 it is also necessary to change some kinetic law. For model 151 they report a "missing step", but since the opposite reaction is part of the model, once again this amounts to adding a reverse reaction to an existing one. Finally, for model 44 they describe that the issue is that some kinetic expression does not depend on one of the reactants of the reaction, making it possible for that reactant's concentration to become negative.

Model inconsistencies studied in

For models 93, 94 and 151, which indeed are flagged by the "Negative" warning, our algorithm correctly adds the missing reverse reactions, directly from the kinetic expressions. The models automatically curated this way do not raise any warning at the end.

For model 44, the automatic curation allows us to get rid of a "K not R" warning by transforming the reaction v3 A + Y cytosol×V m3×A 4 ×Y 2 ×Z 4 /((Ka 4 +A 4 )×((Ky 2 +Y 2 )×(Kz 4 +Z 4 )))

-→ A + Z into Z + A + Y -→ 2 × Z + A
with the same kinetics. However, as expected, the "R not K" warning identified by Kaleta et al. remains, the obtained model is still not well-formed. The same happens with model 143 where indeed a "R not K" warning remains after automatic curation, in accordance with the earlier results.

Conclusion

We have described an algorithm for trying to infer a meaningful reaction system from a system of ordinary differential equations. This algorithm is based on a general consistency condition between the kinetic expression and the structure of a reaction in terms of its reactants, products and inhibitors.

We have shown some general properties enjoyed by the influence graph of the Jacobian sign matrix associated to such well-formed reaction systems. These theoretical results militate for distinguishing between catalysts and inhibitors in the modifiers of a reaction, and for using structural analysis methods before fixing parameter values and going to simulations.

We have also evaluated the capability of our reaction inference algorithm to automatically curate the writing of ODE models with reactions by applying it to the ODE models generated from the SBML models of the curated part of biomodels.net. In particular, we have shown that the inference of well-formed reactions from the ODEs, combined with the inference of hidden molecules corresponding to linear invariants, is sufficient to automatically curate the writing of some ODE models of the cell cycle with consistent reactions. On the whole curated part of the biomodels.net repository, we have shown that our automatic curation method significantly improves the writing of the models with reactions by reducing the number of non well-formed reaction systems from 65% to 29%.

Although the primary concern of SBML is to provide a common format for exchanging models and doing simulations, we believe that stronger consistency conditions should be enforced in SBML to perform structural analyses, and that the strict well-formedness conditions presented in this paper should be verified by non reversible reactions.

  [START_REF] Angeli | A petri net approach to persistence analysis in chemical reaction networks[END_REF]. The DIG of Example 2.10 can be depicted by the following graph:

FIG. 1 .

 1 FIG.1. Minimal cascade model for mitotic oscillations. Cyclin is synthesized at a constant rate (vi) and triggers the transformation

Algorithm 3 . 1 .

 31 (Hidden molecule inference) input: ODE system O over variables {x 1 , . . . , x s }, 1. iteratively replace in O any expression of the form -x + y by y -x, 2. for each expression of the form k -x -y in O where k is a numerical constant or a parameter, and x and y are variables, (a) introduce a new variable z with time derivative ż = -ẋ -ẏ, and functional dependency equation z = k -x -y, (b) substitute any occurrence of k -x -y in O by z, (c) substitute any occurrence of k + v -x -y in O for any expression v, by v + z, (d) substitute any occurrence of k -x + w -y in O for any w, by v + z, 3. for each expression k -x appearing in O where k is a constant or a parameter and x a variable, (a) introduce a new variable z with time derivative ż = -ẋ and functional dependency equation z = k -x, (b) substitute any occurrence of k -x in O by z, (c) substitute any occurrence of k + v -x in O for any expression v, by z + v, output: ODE system O over variables {x 1 , . . . , x s } and hidden molecule variables {z 1 , . . . , z k },

Algorithm 3 . 6 .

 36 (Reaction inference) input: ODE system O over variables for molecular concentrations, partial_has_pos_val test 1. rewrite O into additive normal form 2. compute the set T of all terms appearing in O 3. let R := ∅ 4. for each non-decomposable term t in T , (a) let r := ∅ , p := ∅ , m := ∅ (b) for each variable x where t occurs with integer coefficient c in ẋ in O, i. if c < 0 then r(x) := -c, ii. if c > 0 then p(x) := c, (c) for each variable x such that r(x) = 0 and partial_has_pos_val(t, x), i. r(x) := 1, ii. p(x) := p(x) + 1, (d) for each variable x such that partial_has_pos_val(-t, x), i. m(x) := 1, (e) R := R ∪ {r / m t -→ p}, output: reaction system R. Example 3.7. The model of three reactions of Example 2.10 has one invariant: pMPF + MPF is indeed a constant c (the sum of initial values of pM P F and M P F ) since

Example 3 . 13 .

 313 For instance, given the ODE system ẋ = -2kx = -ẏ, our algorithm infers the reaction 2 × x k×x -→ 2 × y, whereas a Mass Action kinetic reaction model for this system is x 2×k×x -→ y.

"

This paper is an extended version of a communication presented at CMSB'12. The algorithms described in this paper are implemented in the open-source software modeling platform Biocham[1,

2] available at http://lifeware.inria.fr/biocham/ release 3.5. The models used in the experiments are available from http://www.biomodels.net/ release 24

http://lifeware.inria.fr/biocham/

Acknowledgement: This work has been supported by the French OSEO Biointelligence and ANR Biotempo projects. We are grateful to the partners of these projects for fruitful discussions. We also acknowledge the comments of the reviewers which were useful to improve the presentation of this article.


Sébastien Kubicki
email: sebastien.kubicki@enib.fr

Sophie Lepreux

Christophe Kolski

Distributed User Interfaces: Collaboration and Usability Distributed UI on Interactive tabletops: issues and context model

Keywords:

The User Interface distribution can also be applied on interactive tabletops which are connected and more or less remote. This distribution raises issues which concern collaboration (how to distribute the UI to collaborate?); besides, concerning the tangible interaction: which role and appearance (tangible or virtual) must have the objects? In this chapter we describe an extended context model in order to take into account both interactions on a single interactive tabletop and interactions which are distributed and collaborative. The model proposed can, from our point of view, be used to make sure that the usability of the interaction is guaranteed. Indeed, it is essential to know the interaction configuration in order to ensure the usability of the system. The model suggested is illustrated in a case study integrating collaboration and UI distribution. A conclusion gives the limits of the article before a presentation of prospects.

INTRODUCTION

In a world in which everything and everybody are connected, it is possible to envisage connecting different platforms in order to carry out remote collaborations [START_REF] Lepreux | From Centralized interactive tabletops to Distributed surfaces: the Tangiget concept[END_REF] [START_REF] Vanderdonckt | Distributed user interfaces: How to distribute user interface elements across users, platforms, and environments[END_REF][5] [START_REF] Tesoriero | Distributed User Interfaces: Applications and Challenges[END_REF] (see also the chapters written by Garrido Navarro and his collegues, Rädle and his collegues, Harboe and his collegues in this book). Our work concentrates on the connection of interactive tabletops on which the users interact using tangible objects. There are few works in existence relating to connected tables. Very often, some complementary devices are introduced, such as in work of Yamashita et al. [START_REF] Yamashita | Supporting fluid tabletop collaboration across distances[END_REF] which adds a videoconferencing system in order to support collaboration. In our work, we propose the use of one or several tables and some tangible objects in order to support interaction and collaboration. In [START_REF] Lepreux | A step towards the distribution of tangible and virtual objects[END_REF], collaboration scenarios were presented according to different configurations (type of source platform and target platform, distribution strategy (master/slave or autonomous entities), UI distribution type (complete/partial), collaboration type (synchronous/asynchronous)). As suggested in Figure 1, in case of inter-connected tabletops, problematics concerning centralized distribution of UI, as well as network of DUI are various; for instance: how to connect such interaction supports? How to duplicate and extract information? etc. [START_REF] Lepreux | A step towards the distribution of tangible and virtual objects[END_REF]. Afterwards, in [START_REF] Lepreux | From Centralized interactive tabletops to Distributed surfaces: the Tangiget concept[END_REF], a concept of tangible objects called Tangigets is defined and illustrated, making it possible to support distant collaboration between a tangible interactive tabletop and other surfaces. Considering that the context could be a means to ensure the system usability [START_REF] Yang | Does Context Matter? The Impact of Use Context on Mobile Internet Adoption[END_REF], we think that it is interesting to take the context into account in order to try to propose usable surfaces.

In this article we propose to use the context model suggested by Kubicki et al. [START_REF] Kubicki | Simulation in Contexts Involving an Interactive Table and Tangible Objects[END_REF][13], which we widened so as to integrate the specificities of interactive tabletops in order then to widen it even further to integrate distribution characteristics. The following section of the article introduces the TangiSense interactive tabletop as well as the architecture which enables it to manage distribution. Then, we present the proposed context model. An application of this model is then illustrated before concluding and proposing research prospects.

FROM TANGISENSE INTERACTIVE TABLETOP TO DISTRIBUTED SURFACES

This section aims to give a brief presentation of the support of our work which is the TangiSense interactive tabletop. This table has the characteristic of not being tactile, unlike the majority of interactive tabletops present on the market or in the scientific literature. It proposes a direct interaction via tangible objects. Table 1 presents eight different interactive tabletops, each one using a different capture technology.

It shows that each capture technology has its own characteristics. That is why the current interactive tabletops combine technologies. For more information, see Kubicki et al. [START_REF] Kubicki | RFID-driven situation awareness on TangiSense, a table interacting with tangible objects[END_REF]. The first part of this section presents the table while the second one presents the software architecture adopted to support the collaboration and the distribution between surfaces.

The TangiSense Interactive tabletop

The TangiSense interactive tabletop is a prototype which uses RFID technology in order to communicate with tangible objects (equipped with RFID tag(s)). One tile of RFID reading which composes the table is shown in Figure 2(a). It was designed by the RFIdées1 company. For the users, the tabletop looks like a traditional table, because it is a similar size and texture (glass) to a table for everyday use (e.g. due to the technology employed, the users can place their hands on the table without interfering with the system). However, the tabletop has communication capacities via LEDs on its surface which make it possible to display texts (in low resolution) or to define zones (Figure 2(b)).

Figure 2. Tiles composing the tabletop in its first versions

It is possible, according to the applications, to use an external video projector in order to provide a finer grained display and to project directly on tabletop. The tabletop detects the RFID tagged objects and reacts according to them. Moreover, RFID tags offer the possibility to track objects, to store data into objects or to superimpose objects. In order to support the remote collaboration, particular objects, named Tangigets, were defined [START_REF] Lepreux | From Centralized interactive tabletops to Distributed surfaces: the Tangiget concept[END_REF]. Details on the technical aspects of the table can be found in [START_REF] Kubicki | New Human-Computer Interactions Using Tangible Objects: Application on a Digital Tabletop with RFID Technology[END_REF]. A photograph of the table in an experimentation situation with several users is shown in Figure 3.

In its last evolution, TangiSense integrates directly a screen instead of LEDs.

Architecture dedicated to distributed interaction

The software aspects were initially defined for a single table, i.e. even if tables were physically connected, no interaction was envisaged with other platforms (the initial architecture is presented in Kubicki et al. [START_REF] Kubicki | New Human-Computer Interactions Using Tangible Objects: Application on a Digital Tabletop with RFID Technology[END_REF]). From now on, the table architecture has evolved in order to be able to connect several tables and even several platforms, and so allow distribution [6] [11] [16]. This architecture Figure 4 is based on a Multi-Agent System (MAS) which will make it possible to integrate the context adaptation rules. The architecture is fully described in [START_REF] Garbay | Normative Multi-agent Approach to Support Collaborative Work in Distributed Tangible Environments[END_REF] and allows us to consider new distributed interactions in various contexts. (case of tables in interaction)

PROPOSITION OF CONTEXT INTERACTION MODEL TO SUPPORT THE DISTRIBUTION

Context-aware computing [START_REF] Schilit | Context-Aware Computing Applications[END_REF] appeared along with mobile platforms in order to adapt the applications to these new more restricted devices. In these contexts the user was alone on the platform but could switch from one platform to another (in general only one at the same time) 2 . Since then, many definitions and evolutions have covered the concept of context awareness. However, we chose to base our work on the Calvary et al. [START_REF] Calvary | Adaptation des Interfaces Homme-Machine à leur contexte d'usage[END_REF] proposal which defines the context as a triplet <User, Platform, Environment >. From this triplet and a set of definitions from the state of the art, we proposed a context model [START_REF] Kubicki | Contribution à la prise en considération du contexte dans la conception de tables interactives sous l'angle de l[END_REF]. This model was enriched to take into account the specificities of the interactive tabletops [START_REF] Kubicki | Simulation in Contexts Involving an Interactive Table and Tangible Objects[END_REF].

From a User point of view (Figure 5), the most important modification in comparison to the other platforms is the cardinality between User and Platform. Indeed, the interactive tabletops make it possible to work with several users around their surface. The Location of the user relating to the platform is very important and influences the platform (display) itself. Indeed, the context will not be the same if the user position is on one side of the table (e.g. East) or on the opposite side (e.g. West). User posture, sitting or standing, is also important. Finally the interaction style is Post-WIMP and corresponds to competences in the use of tactile technology which can be multiTouch (tactile interactions), or in the use of objects (tangible interaction). The users who are collaborating together can have a social link which will influence their collaboration (parents/children, manager/employee, husband/wife, etc). It would also be necessary to specify which of the tasks are collaborative. It may only be a few of them (partial distribution) or all of them (full distribution). Problems appear concerning the distribution of the users distributed on tabletops: how to communicate, synchronize the tasks, collaborate? Some particular objects have to be defined to solve such problems [START_REF] Lepreux | From Centralized interactive tabletops to Distributed surfaces: the Tangiget concept[END_REF].

Figure 5. User-centric context-awareness model integrating the specifications for the interaction with interactive tabletops, collaboration and the distribution

From a Platform point of view (Figure 6), an attribute was added making it possible to know if the platform is multiUser or not. In order to detect inconsistencies relating to the position of the user (criteria Location, attribute position), an attribute making it possible to know the height of the platform was added. One of the interactive tabletop's characteristics is also the capacity to recognize or interact with a set of objects. We distinguish two object types: virtual objects (which are generally video-projected or displayed) and tangible objects (physical objects placed on the table). From a general point of view, the tangible objects are equipped with tags. These tags can be of a different type: bar-codes, RFID, etc., often stuck under objects so as to enable their identification. According to the capture system, it is possible to vary the number of tags stuck under the object. Thus with three tags stuck under an object, it is possible, for example, to (re)form the shape of the object using software or to detect the direction of rotation of the object (e.g. in the case of RFID technology). The RFID Technology makes it possible to store information directly in the object. Thus the object has a memory enabling it to be completely independent from the table. In the case of virtual objects, it is a videoprojector which is mainly used (it is the case with the first version of Tangi-Sense). These platforms also evolve as regards the work surface and depend on the capture system. Finally, certain platforms can also carry an Artificial Intelligence element which can be centralized or distributed (e.g. TangiSense). A difficulty about the distribution is to connect many platforms, particularly when they have heterogenous characteristics: UI adaptation may bring promising solutions [START_REF] Kubicki | Simulation in Contexts Involving an Interactive Table and Tangible Objects[END_REF].

Figure 6. Platform-centric context-awareness model integrating the specifications for the interaction with interactive tabletops, collaboration and the distribution

Concerning the Environment (Figure 7), a criterion which was added relates to the collective classification of the environment: either (1) the users use the common part of the interactive tabletop, the type of workspace is then common; or (2) they use different parts of the tabletop, each one having their own workspace, in this case the type of workspace is individual [START_REF] Kubicki | RFID-driven situation awareness on TangiSense, a table interacting with tangible objects[END_REF].

Figure 7. Environment-centric context-awareness model integrating the specifications for the interaction with interactive tabletops, collaboration and the distribution

In order to take the collaboration into account, a collaborative environment criterion was added to the Environment characteristics. Indeed, this is the environment conditions which will possibly lead the users to collaborate. This collaborative environment can be distinguished through two characteristics which will influence the interactions. Either the users collaborate by using the table (meaning a single table), in this case the collaboration is considered to be on site; or the collaboration is done on several distant platforms (the users use the table and at least one other distant support such as another table and/or another platform); in this case the environment is one of remote collaboration. The method of collaboration must be also added to the environment in order to know if collaboration is synchronous or asynchronous.

The next section aims to introduce a simulation to show our proposition and models in real context situations.

SIMULATION

The scenario involves a child working on an amusing activity on a first interactive tabletop: Ricardo performs one or several school exercises. In particular, he is currently using the "learning and recognition of colors" application.

This application intended to teach the recognition and learning of colors (only red, yellow, blue and green) to children (aged from 2 to 5 according to the level of difficulty). The scenario is based on the French teaching syllabus for nursery schools. We asked a teacher to imagine one or more scenarios using an interactive tabletop and a set of objects without giving any limits or constraints. The teacher proposed a simple application in which the children have to move a set of objects which have "lost their color" into the suitably colored frame (i.e. a "black and white" bee should be placed inside a yellow frame) [START_REF] Kubicki | Evaluation of an interactive table with tangible objects: Application with children in a classroom[END_REF].

Meanwhile, his parents, Renato and Sofia, interact in another room around a second interactive tabletop (used with another goal, such as family records). The human-machine interface is different on both platforms in this first context. The child may at various times request assistance (remote collaboration) from his parents. They must be based on a common support to facilitate collaboration (e.g. direct or indirect advice to be provided on the color of certain vegetables): the humanmachine interface must be distributed on the two tabletops in this second context.

An illustration is provided in Figure 8. Several attributes of the context have changed, leading to the adaptation of the distributed interaction.

CONCLUSION & PROSPECTS

This chapter has presented the TangiSense tabletop and its software architecture. This architecture allows UI distribution between different tabletops and surfaces in general. This distribution is possible by integrating the intelligence of the distribution in the agents of a multi-agent system (developed with JADE [START_REF] Bellifemine | Developing multi-agent systems with a FIPAcompliant agent framework[END_REF]). Then a global context model has been described; it has been extended by taking the specificities of tabletops into account. The context model proposed in this chapter considers UI distribution on surfaces.

Our research perspectives are the following: to propose a set of adaptation mechanisms integrated into the multi-agent systems, based on the context model proposed; to propose explicit relations between the context model and usability criteria; to validate progressively the context model in various experiments, including remote collaboration.

Figure 1 .

 1 Figure 1. Problematics concerning (a) centralized distribution of UI, (b) network of DUI

Figure 3 .

 3 Figure 3. TangiSense, its objects and users during experimentation

Figure 4 .

 4 Figure 4. Functional architecture of TangiSense [6](case of tables in interaction)

Figure 8 . 2 )

 82 Figure 8. Illustration of the users in remote collaboration (context #2) Such situations are under study and development. Different evaluations are also planned

Table 1 . Five interactive tabletops using different capture technologies Capture technology Representative example

 1

	Magnetic	Actuated Workbench [19]	No	Yes	No
	RFID	TangiSense [12]	Yes	Yes	Yes

www.rfidees.fr a) b)

However some works such as those of Grolaux et al.[START_REF] Grolaux | Detach Me, Assemble Me Like You Work[END_REF] are proposed on interface distribution between two types of platform.

Acknowledgments This research work was partially financed by the "Ministère de l'Education Nationale, de la Recherche et de la Technologie", the "région Nord Pas de Calais", the "Centre National de la Recherche Scientifique", the FEDER, CISIT, and especially the "Agence Nationale de la Recherche" (ANR TTT and IMAGIT projects ANR 2010 CORD 01701). The authors thank also Yoann Lebrun, Emmanuel Adam and René Mandiau for their contribution in this research.