Vincent Acary

Hidde De Jong

Bernard Brogliato

Numerical Simulation of Piecewise-Linear Models of Gene Regulatory Networks Using Complementarity Systems

Keywords: Nonsmooth dynamical systems, gene regulatory networks, piecewise-linear models, Filippov solutions, differential inclusions, complementarity systems Systèmes dynamiques non réguliers, réseaux de régulation génique, modèles linéaires par morceaux, solutions de Filippov, inclusions différentielles, systèmes de complémentarité

Gene regulatory networks control the response of living cells to changes in their environment. A class of piecewise-linear (PWL) models, which capture the switch-like interactions between genes by means of step functions, has been found useful for describing the dynamics of gene regulatory networks. The step functions lead to discontinuities in the right-hand side of the differential equations. This has motivated extensions of the PWL models based on differential inclusions and Filippov solutions, whose analysis requires sophisticated numerical tools. We present a method for the numerical analysis of one proposed extension, called Aizermann-Pyatnitskii (AP)extension, by reformulating the PWL models as a mixed complementarity system (MCS). This allows the application of powerful methods developed for this class of nonsmooth dynamical systems, in particular those implemented in the Siconos platform. We also show that under a set of reasonable biological assumptions, putting constraints on the right-hand side of the PWL models, AP-extensions and classical Filippov (F)-extensions are equivalent. This means that the proposed numerical method is valid for a range of different solution concepts. We illustrate the practical interest of our approach through the numerical analysis of three well-known networks developed in the field of synthetic biology.

Introduction

When confronted with changing environmental conditions, living systems have a remarkable capacity to rapidly adapt their functioning. For instance, the response of a bacterial cell to the depletion of an essential nutrient leads to the upregulation and downregulation of the expression of up to several hundreds of genes [START_REF] Hengge-Aronis | The general stress response in Escherichia coli[END_REF]. The genes encode enzymes, transcription regulators, membrane transporters and other macromolecules playing a role in cellular processes. The control of the adjustment of gene expression levels is achieved by so-called gene regulatory networks, consisting of genes, RNAs, proteins, and their mutual regulatory interactions.

In order to understand how a particular network structure brings about observed changes in gene expression, mathematical models in combination with computer simulation are increasingly used, especially in the emerging field of systems biology [START_REF] Alon | An Introduction to Systems Biology: Design Principles of Biological Circuits[END_REF]. The modeling of gene regulatory networks also plays a prominent role in synthetic biology, which aims at designing a network structure capable of producing a desired gene expression pattern, for instance a robust oscillation or an externally controlled switch [START_REF] Khalil | Synthetic biology: applications come of age[END_REF]. The networks may be constructed de novo or obtained by rewiring a natural regulatory network.

A variety of formalisms are available for modeling gene regulatory networks [START_REF] Bolouri | Computational Modeling of Gene Regulatory Networks -A Primer[END_REF][START_REF] Karlebach | Modelling and analysis of gene regulatory networks[END_REF]. For many purposes, approximate models based on simplifications of classical kinetic models have been proven useful [START_REF] Chaouiya | Logical modelling of gene regulatory networks with GINsim[END_REF][START_REF] Jong | Strategies for dealing with incomplete information in the modeling of molecular interaction networks[END_REF][START_REF] Glass | Logical and symbolic analysis of robust biological dynamics[END_REF]. First, the approximate models are easier to calibrate against experimental data, due to the fact that they reduce the number of parameters and the complexity of the rate equations. This may help relieve what is currently a bottleneck for modeling in systems biology, namely obtaining reliable estimates of parameter values. Second, the simplified mathematical form of the models makes them easier to analyze. Among other things, this makes it possible to single out the precise role of specific subnetworks [START_REF] Abou-Jaoude | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF][START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF] and to analyze the feasibility of control schemes [START_REF] Edwards | Control design for sustained oscillation in a two gene regulatory network[END_REF].

In this paper we look at one particular class of approximate models of gene regulatory networks, so-called piecewise-linear (PWL) models [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF]. The PWL models are systems of coupled differential equations in which the variables denote concentrations of gene products, typically proteins. The rate of change of a concentration at a particular time-point may be regulated by other proteins through direct or indirect interactions. The PWL models capture these regulatory effects by means of step functions that change their value in a switch-like manner at threshold concentrations of the regulatory proteins. The step functions are approximations of the sigmoidal response functions often found in gene regulation.

PWL models with step functions have favorable mathematical properties, which allows for the analysis of steady states, limit cycles, and their stability [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: The case of a negative feedback loop[END_REF][START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF]. The use of step functions, however, leads to discontinuities in the right-hand side of the differential equations, due to the abrupt changes of the value of a step function at its threshold. These discontinuities are sometimes ignored, which is potentially dangerous as it may cause steady states and other important dynamical properties of the system to be missed. In order to deal with the discontinuities, several authors have proposed the use of differential inclusions and Filippov solutions [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF][START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF]. These proposals to extend PWL models to differential inclusions differ in subtle but nontrivial ways, giving rise to systems with nonequivalent dynamics [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF].

Currently, only few computational tools are available to support the analysis of the differential inclusions obtained from Filippov extensions of PWL models. Genetic Network Analyzer (GNA) provides a qualitative analysis of PWL models of gene regulatory networks (e.g., [START_REF] Batt | Genetic Network Analyzer: a tool for the qualitative modeling and simulation of bacterial regulatory networks[END_REF][START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF]). However, the analysis is based on hyperrectangular overapproximations of the differential inclusions proposed in [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF], and it is currently not clear to which extent this introduces artifacts in the analysis. Moreover, the predictions obtained from this analysis are purely qualitative, describing possible transitions between state-space regions rather than giving numerical solutions. Alterna-tively, an algorithm based on the use of steep sigmoidal response functions in combination with singular perturbation theory has been presented [START_REF] Ironi | A computational framework for qualitative simulation of nonlinear dynamical models of gene-regulatory networks[END_REF].

The aim of this paper is to propose a theoretically sound and practically useful method for the numerical simulation of gene regulatory networks described by PWL models. We notably show that the Aizerman & Pyatnitskii extension (see [37, Definition c), page 55] or [START_REF] Aizerman | Foundation of a theory of discontinuous systems. 1 & 2[END_REF]) of PWL models can be reformulated in the framework of complementarity systems or differential variational inequalities [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF][START_REF] Pang | Differential variational inequalities[END_REF]. The Aizerman & Pyatnitskii extension has been introduced in the context of PWL models of gene regulatory networks in [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF], where it is shown that it leads to a more restrictive extension than the standard Filippov extension. The reformulation as a complementarity system allows us to employ the rich store of numerical methods available for these and other classes of discontinuous systems [START_REF] Acary | Time-stepping numerical simulation of switched circuits with the nonsmooth dynamical systems approach[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. Moreover, we show that under two reasonable biological assumptions, posing constraints on the admissible network structures, the different extensions of PWL models that have been proposed, as well as the hyperrectangular overapproximation in [START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF], are equivalent. This means that the numerical simulation approach developed in this paper is valid for a range of different solution concepts for PWL models of gene regulatory networks.

We illustrate the interest of our numerical simulation approach by means of the analysis of three synthetic networks published in the literature: the repressilator [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF], an oscillator with positive feedback [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF], and the IRMA network [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF]. We develop PWL models of these networks, either from scratch or by adapting existing ODE models, and numerically simulate the dynamic response of these networks to external stimuli. The simulations are shown to reproduce known qualitative features of these networks, notably the capability to generate (damped) oscillations for the first two networks, and a switch-on/switch-off response after a change in growth medium for the third. We believe these examples demonstrate that the numerical simulation approach developed in this paper provides a useful extension of the toolbox of modelers of gene regulatory networks.

PWL models of gene regulatory networks 2.1 Definition of PWL models

The dynamics of genetic regulatory networks can be described by piecewise-linear (PWL) differential equation models using step functions to account for regulatory interactions [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF][START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF][START_REF] Mestl | A mathematical framework for describing and analysing gene regulatory networks[END_REF]. In this section we briefly summarize the PWL modeling framework.

We denote by x = (x 1 , . . . , x n ) T ∈ Ω a vector of cellular protein or RNA concentrations, where Ω ⊂ R n + is a bounded n-dimensional hyperrectangular subspace of R n + . For each concentration variable x i , i ∈ {1, . . . , n}, we distinguish a set of constant, strictly positive threshold concentrations {θ 1 i , . . . , θ pi i }, p i > 0. At its thresholds a protein may affect the expression of genes encoding other proteins or the expression of its own gene. We call Θ = i∈{1,...,n},k∈{1,...,pi} {x ∈ Ω | x i = θ k i } the subspace of Ω defined by the threshold hyperplanes.

Definition 1 (PWL model).

A PWL model of a gene regulatory network is defined by a set of coupled differential equations

ẋi = f i (x) = -γ i x i + b i (x) = -γ i x i + l∈Li κ l i b l i (x), i ∈ {1, . . . , n}, (1)
where κ l i and γ i are positive synthesis and degradation constants, respectively, L i ⊂ N are sets of indices of regulation terms, and b l i : Ω \ Θ → {0, 1} are so-called regulation functions.

Inria

Intuitively, (1) defines the rate of change of each concentration x i as the difference of the rate of synthesis (the first term in the right-hand side) and the rate of degradation (the second term). The synthesis term depends on the concentrations of regulatory proteins through the regulation functions, which account for the interactions between the genes in the network. Degradation is described by a first-order term including contributions of growth dilution and protein degradation. While this is sufficient for the examples treated in this paper, the degradation term in (1) can be easily extended to include proteolytic regulators.

Each regulation function b l i (•) is defined in terms of step functions

s + (x j , θ k j ) = 1 if x j > θ k j 0 if x j < θ k j and s -(x j , θ k j ) = 0 if x j > θ k j 1 if x j < θ k j , (2)
where x j is a concentration variable, j ∈ {1, . . . n}, and θ k j a threshold for x j , k ∈ {1, . . . p i }. Notice that s -(x j , θ k j ) = 1 -s + (x j , θ k j ). The step functions capture the switch-like character of gene regulation by transcription factors and other proteins. The regulation functions are algebraic equivalents of discrete Boolean functions expressing the combinatorial logic of gene regulation [START_REF] Thomas | Biological Feedback[END_REF].

For future use, we introduce the following generalizations of the step functions s + (x j , θ k j ) and s -(x j , θ k j ), which consists of extending the function at the discontinuities:

S + (x j , θ k j ) =    1 x j > θ k j [0,1] x j = θ k j 0 x j < θ k j and S -(x j , θ k j ) =    0 x j > θ k j [0,1] x j = θ k j 1 x j < θ k j . (3)
Figure 1 gives an example of a PWL model of a simple two-gene network. The network consists of two genes each of which encodes a protein that inhibits the expression of its own gene and activates the expression of the other gene. Activation and inhibition are assumed to occur at different thresholds. In order for a gene to be expressed, its activator needs to be present (above its threshold) and its inhibitor absent (below its threshold). While this example is used for illustrative purposes in this and the next section, it should be noted that it includes many aspects of actual regulatory networks: auto-regulation, cross-regulation, and combinatorial regulation of gene expression (see Section 6). Several PWL models of actual regulatory networks are available in the literature, see [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF] and references therein.

Protein 1 gene 1 Protein 2 gene 2 (a) ẋ1 = -γ1 x1 + κ1 s + (x2, θ 1 2 ) s -(x1, θ 2 1 ) ẋ2 = -γ2 x2 + κ2 s + (x1, θ 1 1 ) s -(x2, θ 2 2 ) (4) (b)

Regulation functions in PWL models

Which regulation functions b l i (•) entering the PWL models are admissible? In order to answer this question, we first develop how the regulation functions relate to Boolean functions describing the combinatorial control of gene regulation.

Recall that each variable x j has p j thresholds {θ 1 j , . . . , θ pj j }. The step functions can be associated with Boolean variables X k j such that

X k j (x) = (x j > θ k j ) = s + (x j , θ k j ) Xk j (x) = (x j < θ k j ) = s -(x j , θ k j ), (5)
where X denotes the complemented variable of X.

Let us give some basic definitions of Boolean algebra [START_REF] Kohavi | Switching and Finite Automata Theory[END_REF]Chapter 3]. A literal denoted by Y j is defined either as the Boolean variable Y j or its negation Ȳj . Given a set of n Boolean variables Y 1 , . . . , Y n , a minterm m is defined as a conjunction of exactly n literals in which each Y j , j ∈ {1 . . . n} appears once (each variable appears once in either its complemented or uncomplemented form). That is,

m = n j=1 Y j . (6)
For the set of variables X k j , j ∈ {1, . . . , n}, k ∈ {1, . . . , p j }, we have 2 p minterms, with p = j∈{1,...,n} p j , and these are denoted as

m α (x) = n j=1 pj k=1 X k j (x), α ∈ {0, . . . , 2 p -1}. (7)
The subscript α corresponds to the decimal equivalent of the conventional binary encoding of literals (1 for X k j and 0 for Xk j ). Table 1 gives an example of the minterms for the four variables in the example of Figure 1, with their corresponding binary encoding and decimal equivalent.

Regulation functions can now be defined as the minterm expression of a Boolean function [START_REF] Plahte | A methodological basis for description and analysis of systems with complex switch-like interactions[END_REF]. More precisely, they can be rewritten in minterm disjunctive normal form (DNF):

b l i (x) = 2 p -1 α=0 c l i,α m α (x), (8)
with c l i,α ∈ {0, 1}. By means of (8), the PWL models can be written as

ẋi = f i (x) = -γ i x i + l∈Li κ l i 2 p -1 α=0 c l i,α m α (x), i ∈ {1, . . . , n}. (9)
As an example, consider the DNF of the regulation functions in the PWL model of Figure 1. Since L 1 = L 2 = 1, we omit the superscript l for κ l i and c l i,α . As each variable has two thresholds, the canonical representation is composed of the 16 minterms shown in Table 1. Notice that for each state variable x 1 , x 2 , only 4 of the coefficients c i,α equal 1. We thus obtain the following PWL model:

ẋ1 = -γ 1 x 1 + κ 1 15 α=0 c 1,α m α (x) ẋ2 = -γ 2 x 2 + κ 2 15 α=0 c 2,α m α (x), (10)

Inria

Binary encoding decimal equivalent minterm constant X 1

Table 1: Minterms and their binary encoding for the PWL model in Figure 1.

with      c 1,2 = c 1,3 = c 1,10 = c 1,11 = 1 c 2,8 = c 2,10 = c 2,12 = c 2,14 = 1 c 1,α = c 2,α = 0 otherwise. (11)
Writing out the right-hand side of (10) yields

ẋ1 = -γ 1 x 1 + κ 1 (m 2 (x) + m 3 (x) + m 10 (x) + m 11 (x)) ẋ2 = -γ 2 x 2 + κ 2 (m 8 (x) + m 10 (x) + m 12 (x) + m 14 (x)), (12)
which can be easily simplified to (4) by noting that

s + (x 2 , θ 2 2 ) = 1 -s -(x 2 , θ 2 2 ) and s + (x 1 , θ 1 1 ) = 1 -s -(x 1 , θ 1
1 ). More generally, minimization procedures from the field of circuit design and finite automata theory [START_REF] Kohavi | Switching and Finite Automata Theory[END_REF]Chapter 4] can be used to reduce the canonical representation to polynomials that are simpler to handle in practice.

The definition of regulation functions in DNF motivates the following modeling assumption [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF].

Assumption 2. The regulation functions b l i (•) are multiaffine functions, that is, they are affine with respect to each s + (x j , θ k j ), for j ∈ {1, . . . n} and k ∈ {1, . . . , p j }. This assumption can be shown to be generic for all regulation functions corresponding to Boolean functions written in DNF. The proposition directly follows from the observation that the regulation functions in (9) are sums of minterms, each of which is affine with respect to the step functions s + (x j , θ k j ), bearing in mind that s -(x j , θ k j ) = 1 -s + (x j , θ k j ).

A second assumption requires that, when two genes have a common regulation, the latter does not act upon the two genes at the same threshold. Assumption 4. Every step function s + (x j , θ k j ), with j ∈ {1, . . . n} and k ∈ {1, . . . , p j }, occurs in at most one b i (•), i ∈ {1, . . . , n}. As a consequence, for a given j, k, every vector [∂b i (x)/∂s + jk ] i∈{1,...,n} , with s + jk ≡ s + (x j , θ k j ), has at most one non-zero element. Assumption 4 is a rather weak modeling assumption, in the sense that there is usually no compelling biological reason for two genes to be regulated at exactly the same threshold. The notable exception is the case of bacterial genes co-transcribed from the same promoter, that is, genes that are included in the same operon. Notice that all models of biological networks presented in Section 6 satisfy Assumption 4.

The interest of the above restrictions on regulation functions, and thus on the right-hand sides of the PWL models, is that together they entail the equivalency of different solution concepts. This will be shown in the next section.

3 Solutions of PWL models

Filippov extensions of PWL models

The use of step functions s ± (x j , θ k j ) in [START_REF] Aizerman | Foundation of a theory of discontinuous systems. 1 & 2[END_REF] gives rise to mathematical complications, because the step functions are undefined and discontinuous at

x j = θ k j . Therefore, f (•) = (f 1 (•), . . . , f n (•))
T is undefined and may be discontinuous on the threshold hyperplanes Θ. In order to deal with this problem, we can follow an approach originally proposed by Filippov [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF] and widely used in control theory. It consists in extending the differential equation ẋ = f (x), x ∈ Ω \ Θ, to a differential inclusion. Following the book of Filippov, Machina and Ponosov [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF] review several different ways in which this can be done. Below we discuss the two main alternatives, which we call F-and AP-extensions, respectively, and we give precise definitions of the corresponding solutions of the PWL models.

Definition 5 (F-extension of PWL models). The F-extension of the PWL model ( 1) is defined by the differential inclusion

ẋ ∈ F(x), with F(x) = co { lim y→x, y ∈Θ f (y)} , x ∈ Ω, ( 13
)
where co(P ) denotes the closed convex hull of the set P , and {lim y→x, y ∈Θ f (y)} the set of all limit values of f (y), for y ∈ Θ and y → x.

This definition corresponds to the classical Filippov approach [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF], as applied in the context of gene regulatory network modeling in [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF].

Formally, we define a F-PWL system Σ as the triple (Ω, Θ, F), that is, the set-valued function F(•) given by [START_REF] Batt | Genetic Network Analyzer: a tool for the qualitative modeling and simulation of bacterial regulatory networks[END_REF], defined on the n-dimensional state space Ω, with Θ the union of the threshold hyperplanes [START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF]. Definition 6. A solution of an F-PWL system Σ on a time interval I is a solution of the differential inclusion (13) on I, that is, an absolutely-continuous vector-valued function ξ(•) such that ξ(t) ∈ F(ξ(t)) almost everywhere on I.

For all x 0 ∈ Ω and τ ∈ R + ∪ {∞}, Ξ Σ (x 0 , τ ) denotes the set of solutions ξ(t) of the PWL system Σ, for the initial condition ξ(0) = x 0 , and t ∈ [0, τ ]. In particular, notice that the derivative of ξ(•) may not exist, and therefore ξ(t) ∈ F(ξ(t)) may not hold, if ξ reaches or leaves Inria Θ at t. The existence of at least one solution ξ on some time interval [0, τ ], τ > 0, with initial condition ξ(0) = x 0 is guaranteed for all x 0 ∈ Ω [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF]. However, in general there is not a unique solution.

As stated above, the book of Filippov proposes other extensions of the PWL systems, which do not define the right-hand side of the inclusion as the limit values of the function f (y), like in [START_REF] Batt | Genetic Network Analyzer: a tool for the qualitative modeling and simulation of bacterial regulatory networks[END_REF]. A common definition of the right-hand side, following [37, Definition c), page 55], is attributed to [START_REF] Aizerman | Foundation of a theory of discontinuous systems. 1 & 2[END_REF] for systems of the form

ẋ = f (x, u), x ∈ R n , u ∈ R p , (14)
where the function f : R n+p → R n is continuous in the set of arguments and u(x) : R n → R p is discontinuous. Systems ( 14) are often encountered in control theory, especially in variable structure systems and sliding mode control [START_REF] Orlov | Discontinuous Systems[END_REF][START_REF] Utkin | Variable structure systems with sliding modes: A survey[END_REF][START_REF] Utkin | Sliding Modes in Control Optimization[END_REF]. At the point of discontinuity, for i = j, the arguments u i and u j of f are supposed to vary independently on the sets U i (x) and U j (x), usually assumed to be closed convex sets in R. The right-hand side of the differential inclusion in the sense of Aizerman and Pyatnitskii is then defined as

G(x) = y | y = f (x, u), u i ∈ U i (x), i ∈ {1 . . . p} (15)
The set-valued vector field G(x) is generically non-convex and therefore often replaced by its convexification [37, page 56] and [START_REF] Utkin | Sliding Modes in Control Optimization[END_REF]. Furthermore, if the arguments u i , i ∈ {1, . . . , p}, are discontinuous on surfaces S i , such that the surfaces are different and the normal vectors at the points of intersection are not linearly dependent, then

H(x) = co (G(x)). A standard result states that if f (•) in (14) is linear in u, then G(x) = H(x) if all U i , i ∈ {1, . . . , m}, are convex
F(x) = G(x) = H(x).
In the context of the modeling of gene regulatory networks, Machina and Ponosov [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF] apply the alternative extension G of the PWL models, using the generalized step functions S + (x j , θ j ) and S -(x j , θ j ) (Section 2.1) instead of the set-valued functions U i . They argue that this extension gives results that are closer to those obtained with gene regulatory network models using sigmoidal functions rather than step functions, in the limit of infinitely-steep sigmoids [START_REF] Plahte | Analysis and generic properties of gene regulatory networks with graded response functions[END_REF]. Moreover, as we will show below, this definition is more convenient for numerical simulation purposes. We call the resulting differential inclusion the Aizerman-Pyatnitskii (AP)-extension of PWL models.

Remark 7. In the seminal book of Filippov [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF], G(•) is denoted as F 1 (•) and H(•) as F 2 (•). In order to avoid confusion with the components of F(•), we choose the alternative notation proposed here.

Let σ = (σ 1 1 , . . . , σ p1 1 , . . . , σ 1 n , . . . , σ pn n ) T ∈ [0, 1] p . Moreover, define the function g : R p → R n by g i (σ) = l∈Ln κ l i bl i (σ), j ∈ {1, . . . , n}, (16)
where bl i (•) are obtained from b l i (•) by replacing every occurrence of s + (x j , θ k j ) and s -(x j , θ k j ) by σ k j and 1 -σ k j , respectively, for all j ∈ {1, . . . , n} and k ∈ {1, . . . , p j }. Definition 8 (AP-extension of PWL models). The AP-extension of a PWL model ( 1) is defined by the following differential inclusion

ẋ ∈    G 1 (x)
. . .

G n (x)    =         -γ 1 x 1 + g 1 (σ) . . . -γ n x n + g n (σ)    σ k j ∈ S + (x j , θ k j ), j ∈ {1, . . . , n}, k ∈ {1, . . . , p j }      . ( 17
)

RR n°8207

In line with this definition, we obtain an AP-PWL system P given by the triple (Ω, Θ, G). The solutions of this system are defined as follows.

Definition 9. A solution of an AP-PWL system P on a time interval I is a solution of the differential inclusion [START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF] on I, that is, an absolutely-continuous vector-valued function ξ(•) such that ξ(t) ∈ G(ξ(t)) almost everywhere on I.

The set of solutions of the AP-PWL system P is denoted by Ξ P (x 0 , τ ), for ξ(0) = x 0 , and t ∈ [0, τ ]. Since G(•) may not be convex contrary to F(•), the existence of solutions of this system for every x 0 ∈ Ω cannot be guaranteed in general.

Notice that the definition of the AP-extension requires that all occurrences of a step function s + (x j , θ k j ) in ( 17) are replaced by the same value σ k j ∈ S + (x j , θ k j ), and all occurrences of s -(x j , θ k j ) by the same value 1 -σ k j . In other words, all occurrences of a step function s + (x j , θ k j ) (respectively s -(x j , θ k j )) are replaced by a selection σ k j (respectively 1 -σ k j ). From a modelling point of view this makes sense, as the different occurrences in b(•) of the step function s ± (x j , θ k j ) usually correspond to a single underlying biophysical process. However, from a mathematical point of view one could imagine to relax the definition by allowing different values σ k j ∈ S + (x j , θ k j ) for different occurrences of s + (x j , θ k j ), and by decoupling the values for positive and negative step functions. This corresponds to replacing [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF] by the alternative definition Ĝ

(x) = f (x, U (x)), where U (x) = [U i (x)] T {i∈1...p} .
In A we show that, in general, this notation is ambiguous and leads to an extension of the PWL systems that is not equivalent with Definition 8. This ambiguity is discussed in [70, Section 1.3], where two inputs u i and u j in a controlled system of the form ( 14) are subjected to discontinuities over the same surface S i = S j .

Relations between different Filippov extensions of PWL models

The natural question to ask is how the solutions of the F-PWL system given by the differential inclusions [START_REF] Batt | Genetic Network Analyzer: a tool for the qualitative modeling and simulation of bacterial regulatory networks[END_REF] relate to the solutions of the AP-PWL system given by [START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF]. This means that we have to compare F(•) and G(•). Proposition 10. Under Assumption 2, F(x) = co(G(x)) for all x ∈ Ω. This result has been proven in [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF]. Since in general G(x) is not a convex set (G(x) = co(G(x))), it follows from the proposition that the two solution concepts may give different results. We illustrate this by means of the following example from [START_REF] Machina | Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks[END_REF]:

ẋ1 = -γ 1 x 1 + κ 1 [s + (x 1 , θ 1 ) + s + (x 2 , θ 2 ) -2 s + (x 1 , θ 1 ) s + (x 2 , θ 2 )] ẋ2 = -γ 2 x 2 + κ 2 [1 -s + (x 1 , θ 1 ) s + (x 2 , θ 2 )]. (18)
Consider this PWL system at the intersection of the two thresholds, that is, at x = (θ 1 , θ 2 ) T . According to Definition 5, we have

F(x) = -γ 1 x 1 -γ 2 x 2 + co 0 0 , 0 κ 2 , κ 1 κ 2 , ( 19
)
that is, the convex combination of the vector fields in the four regions having (θ 1 , θ 2 ) T in their boundary, evaluated at this point. Notice that the vector fields in the regions

{x 1 > θ 1 , x 2 < θ 2 } and {x 1 < θ 1 , x 2 > θ 2 } are the same ([-γ 1 x 1 + k 1 , -γ 2 x 2 + k 2 ] T
), which explains that only three (instead of four) vectors appear in [START_REF] Brogliato | Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings[END_REF]. Figure 2(a) shows the convex envelope of the vector fields at the intersection of the thresholds, for the case that κ 1 > γ 1 θ 1 and κ 2 > γ 2 θ 2 .

Inria

The AP-extension of the PWL model ( 18), according to Definition 8, is defined as follows:

G(x) = -γ 1 x 1 + κ 1 [σ 1 + σ 2 -2σ 1 σ 2 ] -γ 2 x 2 + κ 2 [1 -σ 1 σ 2 ] σ 1 , σ 2 ∈ [0, 1] . ( 20
)
As can be seen in as defined in [START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF]. It includes the vector fields in the neighboring regions of x = (θ 1 , θ 2 ) T , which are obtained when σ 1 , σ 2 take their extreme values 0 or 1. The plots are obtained for κ 1 > γ 1 θ 1 and κ 2 > γ 2 θ 2 .

The example shows that the different Filippov extensions are not equivalent, but the question can be posed if differences occur when analyzing biologically relevant network structures. For example, the network represented by the PWL model ( 18) consists of two regulators that jointly regulate their own genes, at exactly the same threshold concentrations, according to an XOR switch in the case of x 1 and a NAND switch in the case of x 2 .1 A priori this is a rather unlikely configuration to occur in real biological networks.

The following result shows that under Assumptions 2 and 4, which are usually not restrictive in practice, the different solution concepts proposed in the previous section coincide.

Proposition 11. Under Assumptions 2 and 4, F(x) = G(x) for all x ∈ Ω.

Proof. We will prove the proposition by showing that G(x) = co(G(x)), so that the equality of F(x) and G(x) directly follows from Prop. 10.

Following Definition 8, G i (x) can be written in terms of g i (σ), with σ ∈ [0, 1] p and i ∈ {1, . . . , n}. In what follows, it will be convenient to rewrite g i (σ) = g i (ρ i ), where ρ i is the vector including all elements in σ appearing in b i (•), that is, those elements in σ occurring in the right-hand side of [START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF]. Due to Assumption 4, ρ i and ρ j , i = j, have an empty intersection.

We now show by mathematical induction on the dimension k of ρ i , k ≥ 1, that g i (ρ i ) is a closed interval. For k = 1, it follows from Assumption 2 that g i (ρ 1 i ) can be rewritten as

a 0 i + a 1 i ρ 1 i , where a 0 i , a 1 i ∈ R. Since ρ 1 i ∈ [0, 1], ρ 1 i = 0, or ρ 1 i = 1, g i (ρ 1 i ) is a closed interval. Suppose that for some fixed k > 1, g i (ρ 1 i , . . . , ρ k i ) is a closed interval (induction hypothesis). Then, by Assumption 2, g i (ρ 1 i , . . . , ρ k+1 i
) can be rewritten as

g 0 i (ρ 1 i , . . . , ρ k i ) + g 1 i (ρ 1 i , . . . , ρ k i ) ρ k+1 i , for some g 0 i (•) and g 1 i (•).
From the induction hypothesis it follows that the expressions g 0 i (ρ 1 i , . . . , ρ k i ) and

g 1 i (ρ 1 i , . . . , ρ k i ) denote closed intervals, so with ρ k+1 i ∈ [0, 1], ρ k+1 i = 0, or ρ k+1 i = 1, we conclude that g i (ρ 1 i , . . . , ρ k+1 i
) is a closed interval. This proves the assertion. Due to Assumption 4 and the above, all g i (ρ i ) and therefore G i (x) can be independently varied, i ∈ {1, . . . , n}. In other words, G(x) can be written as G 1 (x)×. . .×G n (x) and is a closed hyperrectangle. As a consequence, G(x) is a closed convex set and G(x) = co(G(x)) = F(x) using Proposition 10.

The PWL model of Figure 1 satisfies Assumptions 2 and 4, so that for this example the two Filippov extensions are equivalent. This can be illustrated for a point x on the segment {θ 1 1 < x 1 < θ 2 1 , x 2 = θ 2 2 } (see Figure 3). A straightforward computation yields:

ẋ ∈ F(x) = -γ 1 x 1 + κ 1 -γ 2 x 2 + κ 2 σ 2 σ 2 ∈ [0, 1] . (21)
By means of Definition 8, the same set is obtained for G(x).

= (x 1 , x 2 ) on the segment {θ 1 1 < x 1 < θ 2 1 , x 2 = θ 2 2 }.
In this example, F(x) and G(x) coincide.

As an aside, we notice that the following corollary holds, which makes explicit an intermediate result in the proof of Prop. 11.

Corollary 12. Under Assumptions 2 and 4, F(x) and G(x) are hyperrectangular sets.

The qualitative analysis developed in [START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF][START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF], and implemented in the computer tool GNA, exploits the mathematical properties of hyperrectangular overapproximations of the convex sets Inria used in the classical Filippov approach. In particular, they facilitate the construction of discrete abstractions of the continuous dynamics, in the form of state transition graphs, in the absence of precise quantitative information on the model parameters [START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF][START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF]. The above corollary shows that under reasonable and verifiable modeling assumptions, the results obtained are the same as those obtained with the F-and AP-extensions of Definitions 5 and 8, respectively. In other words, under the assumptions of Prop. 11, the hyperrectangular sets are not overapproximations but exact.

While a qualitative analysis is appropriate for certain problems, the absence of precise quantitative predictions is not desirable in others, such as the analysis of a limit cycle [START_REF] Abou-Jaoude | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF] or the design of a controller for a synthetic network [START_REF] Edwards | Control design for sustained oscillation in a two gene regulatory network[END_REF]. The quantitative study of PWL models of gene regulatory networks is hindered by the fact that, as of now, no general and efficient tool for the numerical simulation of Filippov extensions is available. In the remainder of this paper we show how tools developed for the simulation of nonsmooth mechanical, electrical and control systems [START_REF] Acary | Time-stepping numerical simulation of switched circuits with the nonsmooth dynamical systems approach[END_REF][START_REF] Acary | Nonsmooth Modeling and Simulation for Switched Circuits[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Acary | Implicit Euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF] can be adapted for this purpose.

Numerical methods for AP-extensions of PWL systems

In this section, we propose numerical methods for performing the time-integration of the APextension in Definition 8 and for computing its equilibrium points. General results of convergence (and existence) of solutions are beyond the scope of this paper, but we prove under the assumption that g(•) is continuous (satisfied under Assumption 2) that the discrete one-step problem is always solvable. In practice, we can therefore always numerically compute a solution of the time-discretized problem. In other words, this enables the computation of a selection in the set-valued AP-extension G(•) for the time-discrete system.

Reformulation of PWL models as Mixed Complementarity Systems (MCS)

The main features of the proposed numerical time-integration method are the following:

1. A reformulation of the set-valued relation

σ ∈ S + (x, θ) (22)
in terms of well-known concepts from Convex Analysis and Optimization Theory (inclusion into normal cones, Complementarity Problems (CP) and finite-dimensional Variational Inequalities (VI), see [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF]). One of the interests for introducing such concepts is the extensive amount of work that has been done for their mathematical analysis and their numerical treatment (see [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Rockafellar | Variational Analysis[END_REF] for an overview).

2. An implicit event-capturing time-stepping scheme, mainly based on the backward Euler scheme, which allows to deal with the switch-like behaviour of the flow without resorting to an accurate detection of events. Furthermore, when a sliding motion occurs, such a scheme avoids the "chattering" effect when a attractive surface is reached [START_REF] Acary | Implicit Euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF].

3. The use of efficient numerical solvers for the one-step problem resulting from the timediscretization of the CP/VI formulation. At each time-step, we have to solve a CP (or equivalently a VI) for which numerous efficient solvers exist [START_REF] Billups | Algorithms for Complementarity Problems and Generalized Equations[END_REF][START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF][START_REF] Dirkse | The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems[END_REF][START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Ferris | Interfaces to PATH 3.0: Design, implementation and usage[END_REF]. Efficient enumerative algorithms may also be used to list all possible solutions. When we are concerned with following a single solution trajectory, standard solvers for CP/VI succeed in doing this in an efficient way.

Let us start with the equivalent formulation of ( 22) as an inclusion in the normal cone to the interval [0, 1]. We recall the definition of the normal cone to a convex set C at a point σ ∈ C:

N C (σ) = {x | x T (σ -σ) ≤ 0 for all σ ∈ C}. ( 23
)
By noting that

N [0,1] (σ) =    R -σ = 0 0, σ ∈ [0, 1] R + σ = 1, (24)
the relation

σ ∈ S + (x, θ) (25)
can equivalently be reformulated as the following inclusion

(x -θ) ∈ N [0,1] (σ). ( 26
)
In turn, the relation ( 26) is equivalent to the complementarity conditions

0 ≤ 1 -σ ⊥ (x -θ) + ≥ 0 0 ≤ σ ⊥ (x -θ) -≥ 0, ( 27
)
where the symbol x ⊥ y means x T y = 0 and y + , y -stand for the positive and negative parts of y, respectively. Finally, an equivalent formulation of ( 26) is given by the following VI : find

σ ∈ [0, 1] such that (θ -x) T (σ -σ ) ≥ 0 for all σ ∈ [0, 1]. (28)
Let us define the notion of modes for the solution of the inclusion [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF], or equivalently for the CP/VI formulations. Each solution is associated to a triplet of index sets

I + = {(j, k) ∈ {1 . . . n} × {1 . . . p j } | σ k j = 1, x j -θ k j ≥ 0} I 0 = {(j, k) ∈ {1 . . . n} × {1 . . . p j } | σ k j ∈ (-1, 1), x j -θ k j = 0} I -= {(j, k) ∈ {1 . . . n} × {1 . . . p j } | σ k j = -1, x j -θ k j ≤ 0}, (29)
such that card(I + ∪ I -∪ I 0 ) = p. A mode M α ⊂ Ω × IR p , α ∈ {1 . . . 3 p } is defined by the set of couples (x, σ) such that the triplet of index-sets is equal to a given constant triplet of index sets. Since, in IN p , there are 3 p possibilities to choose these triplets, the total number of modes is 3 p . Let us now define the affine function y : R n → R p such that

y(x) = Cx -θ =              x 1 -θ 1 1 . . . x 1 -θ p1 1 . . . x n -θ n 1 . . . x n -θ pn 1              T ∈ R p , (30)
where C ∈ R p×n , with C ij ∈ {0, 1}, and θ = [θ 1 1 , . . . , θ p1 1 , . . . , θ 1 n , . . . , θ pn n ] T . Using the formulation [START_REF] Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF] and the definition of y(•), the AP-extension of the PWL system [START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF] in Definition 8 can be written compactly as

ẋ = -diag(γ)x + g(σ) y(x) = Cx -θ ∈ N [0,1] p (σ) (31)
Inria where diag(γ) ∈ R n×n is the diagonal matrix composed of γ i , i ∈ {1 . . . n}.

When rewritten in the form (31), the system is a Mixed Complementarity System (MCS). MCSs are an extension of Linear Complementarity Systems (LCS) [START_REF] Heemels | Linear complementarity systems[END_REF] to the nonlinear case with nontrivial bounds on σ. LCS and MCS enter in the more general framework of Differential Variational Inequalities [START_REF] Pang | Differential variational inequalities[END_REF] where the solution set of a variational inequality comes into play in the right-hand side of an ordinary differential equation. For more details on complementarity systems and their relations with other types of nonsmooth dynamical systems, we refer to [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF][START_REF] Georgescu | Switching, relay and complementarity systems: A tutorial on their well-posedness and relationships[END_REF].

Remark 13 (Equilibrium (stationary) points). Finding equilibrium points of the AP-extension of PWL systems is equivalent to solving the following MCP

0 = -diag(γ)x + g(σ) y(x) = Cx -θ ∈ N [0,1] p (σ) ( 32
)
or more compactly

Cdiag(1/(1 + γ))g(σ) -θ ∈ N [0,1] p (σ). ( 33
)
Using the fact that σ → Cdiag(1/(1 + γ))g(σ) -θ is a continuous function and [0, 1] p is a compact convex set of R p , the VI/CP (33) has a nonempty compact set of solutions (direct application of Corollary 2.2.5 in [34, page 148]).

Equilibrium points can then be computed by solving the following problem

     Cdiag(1/(1 + γ))g(σ) -θ ∈ N [0,1] p (σ) x = diag(1/(1 + γ))g(σ) x ∈ R n + . (34)

The general time-discretization framework.

Numerical time-integration of an MCS (31) can be performed by two main families of solvers. The first one, often called event-detecting (or event-driven) time-stepping schemes, performs an accurate location of the time of events. An event corresponds to a change of mode of solutions as defined in the previous section. Between two events, any algorithm for Differential Algebraic Equations (DAEs) can be used for the time integration of the smooth dynamics. In addition to the problem of the drift of constraints when the trajectory is sliding, which can render the detection of events difficult, the main drawback of event-detecting schemes is the inability to deal with infinite accumulations of events in finite time (also termed as Zeno-phenomena). In the context of ordinary differential equations with a discontinuous right-hand side, an eventdriven scheme for a class of piecewise-smooth systems reformulated as complementarity systems is rigorously described in [START_REF] Stewart | A high accuracy method for solving ODEs with discontinuous rigth-hand-side[END_REF]. For an overview of event-driven schemes for nonsmooth systems, we refer to [4, Chapter 7].

The second family of schemes is the class of the event-capturing time-stepping schemes. In this case, no accurate detection of events is performed and the events may occur within the time-step. Although these methods are of low order, they are robust, stable and enjoy some convergence results under specified assumptions (see the survey paper [START_REF] Dontchev | Difference methods for differential inclusions: a survey[END_REF] in the context of general differential inclusions). Moreover, they are able to deal with an infinite number of events in finite time, which is common in practice (see the example in Section 6.2).

For our specific class of inclusions, we are interested in developing a dedicated event-capturing time-stepping scheme which takes benefits from the special structure of the complementarity systems. As we said before, for such systems, the discretization results at each time-step in a CP/VI for which there exists a large amount of efficient algorithms. Therefore, computing a selection of the right-hand side set-valued map of our time-discretization differential inclusion can be done in a very efficient way. For general studies and a survey on time-stepping schemes for complementarity systems, we refer to [START_REF] Acary | Higher order Moreau's sweeping process: mathematical formulation and numerical simulation[END_REF][START_REF] Camlibel | Consistency of a time-stepping method for a class of piecewise-linear networks[END_REF][START_REF] Han | Convergence of time-stepping schemes for passive and extended linear complementarity systems[END_REF][START_REF] Pang | Differential variational inequalities[END_REF] and [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]Chapter 9].

Let us now expose the proposed time discretization of (31) over a time-interval [t k , t k+1 ] of length h:

     x k+1 = x k -h diag(γ)x k+τ + h g(σ k+1 ) y k+1 = Cx k+1 -θ y k+1 ∈ N [0,1] p (σ k+1 ), (35)
with the initial condition x 0 = x(t 0 ). The notation

x k+τ means τ x k+1 + (1 -τ )x k for τ ∈ [0, 1].
The problem ( 35) is called the one-step problem that we have to solve on each interval. Let us define the vector

z k+1 = x k+1 σ k+1 ∈ R n+p , (36)
and the function F : R n+p → R n+p as

F (z k+1 ) = x k+1 -x k + h diag(γ)x k+τ -h g(σ k+1 ) θ -Cx k+1 . ( 37
)
Then the problem ( 35) is equivalent to the following inclusion

-F (z k+1 ) ∈ N R n ×[0,1] p (z k+1 ). (38)
Proposition 14. Let F : R n+p → R n+p be the function defined in [START_REF] Filippov | Differential Equations with Discontinuous Right Hand Sides[END_REF]. Under Assumption 2, the problem of finding

z ∈ R n × [0, 1] p such that -F (z) ∈ N R n ×[0,1] p (z) (39)
has a nonempty and compact solution set.

Proof. The proof is based on the application of Corollary 2.2.5 in [34, page 148] which states that the VI y(σ)

T (σ -σ ) ≥ 0, for all σ ∈ K ⊂ R p , (40)
has a solution if K is compact convex, and y : K → R p is continuous. The inclusion (38) can be restated as

x -x k + h diag(γ)(τ x + (1 -τ )x k ) -h g(σ k+1 ) = 0 Cx -θ ∈ N [0,1] p (σ). (41)
Since the matrix I n + hτ diag(γ) is regular, we can solve the first equation for x. This yields the following definition of y(•):

y(σ) = Cx -θ = C diag(1/(1 + hτ γ)) [(I n -h(1 -τ ) diag(γ))x k + h g(σ)] -θ, (42)
where diag(1/(1+hτ γ)) is the diagonal matrix made of the components 1/(1+hτ γ i ), i ∈ {1 . . . n}.

The inclusion in [START_REF] Glass | Logical and symbolic analysis of robust biological dynamics[END_REF] becomes

y(σ) ∈ N [0,1] p (σ) ( 43
)
and then is equivalent to [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF]. Since g is multiaffine, y is continuous. Choose K = [0, 1] p which is compact convex. Then, the VI (40) has a nonempty and compact solution set for σ. The relation

x = diag(1/(1 + hτ γ)) [(I n -h(1 -τ ) diag(γ))x k + h g(σ)] (44)
allows to build a nonempty and compact set of solutions for z = [x σ] T for [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF].

Inria

Remark 15. The following observations can be made :

• Proposition 14 is crucial for the success of the numerical scheme which is proposed in this paper. It states that we can always find a solution of the discrete problem. It remains to be discussed how to effectively compute one of these solutions.

• Only Assumption 2 is necessary to prove the result. It could even be relaxed to continuous functions g(•) instead of multiaffine.

• To avoid any misunderstanding, we have said nothing about the convergence of the scheme and we cannot straightforwardly extrapolate Proposition 14 for an existence result in continuous time. As explained above, F-extensions of PWL systems ensures existence of solutions.

The AP-extension in general does not guarantee such a result.

At each step, we have to solve the so-called one-step problem which has been proved to possess solutions. In the sequel, we say a few words on the solution procedures in practise. To this end, let us now introduce a fairly standard problem in Mathematical Programming, which is equivalent to the problem [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: The case of a negative feedback loop[END_REF].

Definition 16 (Mixed Complementarity Problem(MCP) [START_REF] Dirkse | The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems[END_REF]). Given a function F : R n+p → R n+p and lower and upper bounds l, u ∈ R n+p ∪{-∞, +∞}, the Mixed Complementarity Problem (MCP) is to find z ∈ R n+p and w, v ∈ R n+p

+ such that        F (z) = w -v l ≤ z ≤ u (u -z) T v = 0 (z -l) T w = 0. ( 45
)
By choosing l ∈ {-∞} n × 0 p and u ∈ {+∞} n × 1 p , solving the inclusion [START_REF] Georgescu | Switching, relay and complementarity systems: A tutorial on their well-posedness and relationships[END_REF] is equivalent to solve an MCP for z k+1 where the vectors w k+1 , v k+1 ∈ R n+p are of the form

(w k+1 ) i = 0, 1 ≤ i ≤ n (y - k+1 ) i-n , n ≤ i ≤ n + m , (v k+1 ) i = 0, 1 ≤ i ≤ n (y + k+1 ) i-n , n ≤ i ≤ n + m . (46)
Numerical algorithms for solving such problems benefit from a long experience (see [START_REF] Billups | Algorithms for Complementarity Problems and Generalized Equations[END_REF][START_REF] Dirkse | The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems[END_REF] for the PATH solver and [START_REF] Billups | A comparison of large scale mixed complementarity problem solvers[END_REF] for a comprehensive comparison of algorithms).

An enumerative procedure based on the Newton-Josephy approach

In most of the applications, the numerical algorithms cited in the previous section are sufficient for computing a solution of the one-step problem. Nevertheless, it can be interesting to be able to enumerate the solutions when they are multiple and isolated. To this end, we can benefit from the fact that y(•) is a linear function of its arguments and g(•) is a multiaffine application for easily applying a Newton-like method. For solving the one-step problem [START_REF] Georgescu | Switching, relay and complementarity systems: A tutorial on their well-posedness and relationships[END_REF] in its MCP formulation, the Newton-Josephy method [START_REF] Josephy | Newton's method for generalized equations[END_REF][START_REF] Rutherford | Miles: A mixed inequality and nonlinear equation solver[END_REF] is used which consists in a linearization of the first line of the time-discretization [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: The case of a negative feedback loop[END_REF]. Then, the MCP is solved by a sequential Mixed Linear Complementarity Problem (MCLP) method using a dedicated enumerative solver for MLCP. Let us introduce the following nonlinear residue

R(x, σ) = (I n + τ hdiag(γ))x -(I n -h(1 -τ )diag(γ))x k -h g(σ). (47)
The solutions x k+1 and σ k+1 of (38) satisfy R(x k+1 , σ k+1 ) = 0. The solutions of the nonlinear system (47) are sought as the limit of the sequence {x ν , σ ν } ν∈N verifying

R L (x ν+1 , σ ν+1 ) = 0 x 0 = x k , σ 0 = σ k ,
where the Newton linearization of R is given by

R L (x ν+1 , σ ν+1 ) = R(x ν , σ ν ) + ∇ x R(x ν , σ ν )(x ν+1 -x ν ) + ∇ σ R(x ν , σ ν )(σ ν+1 -σ ν ).
Let us introduce the so-called iteration matrices

M = ∇ x R(x, σ) = I n + hτ diag(γ) B(σ) = ∇ σ R(x, σ) = ∇ σ g(σ).
At each time-step k, we have to solve the following linearized problem

x ν+1 = M -1 (I n -h(1 -τ )diag(γ))x k + hg(σ ν ) + hB(σ ν )(σ ν+1 -σ α ) . (48)
Inserting into (48) the definition of y(•) given by ( 30), we get the following linear relation between y ν+1 and σ ν+1

y ν+1 = CM -1 (I n -h(1 -τ )diag(γ))x k + hg(σ ν ) + hB(σ ν )(σ ν+1 -σ ν ) -θ.
To summarize, the problem to be solved at each Newton iteration is

y ν+1 = W ν+1 σ ν+1 + q ν+1 y ν+1 ∈ N [0,1] p (σ ν+1 ), (49)
where W ∈ R p×p and q ∈ R p are defined as

W ν+1 = h C M -1 B(σ ν ) q ν+1 = CM -1 [(I n -h(1 -τ )diag(γ))x k + hg(σ ν ) + hB(σ ν )σ ν ] -θ.
The problem ( 49) is an MLCP which can be solved under suitable assumptions by many linear complementarity solvers such as pivoting techniques, interior point techniques and splitting/projection strategies [START_REF] Cottle | The Linear Complementarity Problem[END_REF]. Among these techniques, some efficient enumerative solvers can also be used (see for instance [START_REF] Faiz | An implicit enumeration procedure for the general linear complementarity problem[END_REF][START_REF] Júdice | On the solution of NP-hard linear complementarity problems[END_REF][START_REF] Sherali | Enumeration approach for linear complementarity problems based on a reformulation-linearization technique[END_REF]), if we are interested in enumerating several solutions corresponding to various modes of the switching functions. This approach has therefore an interest from the qualitative point of view since it allows to outline multiple solutions and to choose following one particular solution in the numerical time integration. This feature will be illustrated in the oscillator example of Section 6.2.

Software aspects

The simulation platform Siconos2 developed at INRIA [START_REF] Acary | Nonsmooth Modeling and Simulation for Switched Circuits[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Acary | An introduction to Siconos[END_REF] implements the numerical methods for AP extensions of PWL models described in the previous section. The Siconos platform is an open-source software for the modeling, simulation, analysis and control of nonsmooth dynamical systems. It has been designed and developed with a constant effort to be sufficiently general and modular to be able to deal with applications ranging over several domains, including, Mechanics, Electronics, Control, Systems Theory (see the Siconos website for examples).

In the following, we use the Siconos platform for the simulation of gene regulatory networks. In order to achieve this, the following general formulation of MCS is used

   M ẋ = f (x, t) + g(x, t, λ) y = h(x, t, λ) -y ∈ N [l,u] (λ) (50)
where M is a non-necessarily regular matrix and f (•), g(•) and h(•) are assumed to be userdefined smooth functions. Without any optimization for to the special case of AP-extensions given by ( 31), the computation time for each trajectory is between 15 ms and 1 s on a laptop computer (Apple MacBook Pro 3.60GHz Intel Core 2 Duo). All examples presented in the next section are included in the collection of examples distributed with Siconos.

Illustration of the properties of the numerical methods on a two-gene regulatory network

This section is devoted to the presentation of the properties of the above numerical methods on the genetic regulatory network of two genes in Figure 1, modeled by the two-dimensional system:

ẋ1 = -γ 1 x 1 + κ 1 s + (x 2 , θ 1 2 ) s -(x 1 , θ 2 1 ) ẋ2 = -γ 2 x 2 + κ 2 s + (x 1 , θ 1 1 ) s -(x 2 , θ 2 2 ), (4)
with the following parameters:

θ 1 1 = θ 1 2 = 4, θ 2 1 = θ 2 2 = 8, κ 1 = κ 2 = 40, γ 1 = 4.
5 and γ 2 = 1.5. In particular, we detail the MCS formulation and we propose an analysis in a subset of Ω where all surfaces are attractive. We show that the proposed formulation enables one to show the uniqueness of the solution of the one-step problem and the finite-time stability of one of the equilibrium points.

MCS formulation and stability analysis

We have card(L) = 2, n = 2 and p = 4. Let y(

x) = x 1 -θ 1 1 , x 1 -θ 2 1 , x 2 -θ 1 2 , x 2 -θ 2 2
T and σ ∈ R 4 . We also have

g(σ) = κ 1 σ 1 2 (1 -σ 2 1 ) κ 2 σ 1 1 (1 -σ 2 2 ) , B(σ) = 0 -κ 1 σ 1 2 κ 1 (1 -σ 2 1 ) 0 κ 2 (1 -σ 2 2 ) 0 0 -κ 2 σ 1 1 , M = I 2 + hτ γ 1 0 0 γ 2 , C =     1 0 1 0 0 1 0 1     . (51)
There are three equilibrium points, two of which are attractive ((0, 0) and (θ 2 1 , θ 2 2 )) and a third one repulsive ((θ 1 1 , θ 1 2 )). We focus on the equilibrium point

x 1 = θ 2 1 and x 2 = θ 2 2 by restricting the domain of interest to Ω = Ω ∩ (θ 1 1 , +∞) × (θ 2 1 , +∞).
In Ω, we have σ 1 1 = 1 since x 1 > θ 1 1 , and σ 1 2 = 1 since x 2 > θ 1 2 , so that the original system (31) can be reduced to

ẋ = -diag(γ)x + B σ + κ Cx -θ ∈ N [0,1] 2 (σ), ( 52
) with σ = σ 2 1 σ 2 2 , C = 1 0 0 1 , θ = θ 2 1 θ 2 2 , B = -κ 1 0 0 -κ 2 , κ = κ 1 κ 2 . ( 53
)
The complementarity systems (or equivalently, differential variational inequalities) framework provides not only an efficient framework for simulation, but it can also help us to conclude on the stability, possibly in finite-time, of the equilibrium point (θ 2 1 , θ 2 2 ) in Ω. On one hand, the study of the stability of PWL systems for gene regulatory networks has already been studied in [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF], while on the other hand there is a huge literature dealing with the stability of Filippov systems and differential inclusions, and in particular of the finite-time stability [START_REF] Orlov | Discontinuous Systems[END_REF]. We do not wish to give an exhaustive treatment of the stability of PWL systems by means of monotone differential inclusions, which is beyond the scope of the paper. The following result, however, gives an idea of the large range of applicability of the complementarity systems techniques.

Lemma 17. The equilibrium point θ = (θ 2 1 , θ 2
2 ) T is finite-time stable in Ω. Proof. Since -B is a symmetric definite positive matrix, we choose a symmetric positive definite matrix R as a square root of the inverse of B, i.e. R 2 = -B-1 . Following [START_REF] Brogliato | Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings[END_REF], let us perform the state transformation z = R(x -θ). The solution of the reduced system ( 52) is equivalent to the following inclusion

-ż ∈ T (z), ( 54
)
with the multi-valued operator T defined as

T (z) = R diag(γ)R -1 z+, (diag(γ) θ -κ) + R -1 S + (R -1 z, 0). ( 55
) Since S + (•, 0) is monotone and R = R T , R -1 S + (R -1 •, 0) is also monotone [62, Exercice 12.4].
Since R is a diagonal matrix, we have R diag(γ)R -1 = diag(γ) which obviously is a positive definite matrix. As a consequence, the operator T (•) is a strongly monotone mapping in Ω, i.e., there exists c > 0

(λ -λ) T (x -y) ≥ c x -y 2 , ∀x, y ∈ Ω, ∀λ ∈ T (x), ∀λ ∈ T (y). ( 56
)
Let us compute the set T (0)

T (0) = -diag(γ) θ + κ + BS + ( θ, θ) = [-36, 4] × [-12, 28]. ( 57
)
The existence of an equilibrium point for x = θ is ensured by 0 ∈ T (0). Choosing a Lyapunov function V (z) = 1/2z 2 , we have V (0) = 0, V (z) ≥ 0 for all z and

V (z) = z T ż with -ż ∈ T (z). ( 58
)
Since T (•) is strongly monotone and 0 ∈ T (0), we can conclude that V (z) ≤ c z 2 , c > 0 and that the equilibrium point is exponentially stable. To prove the finite-time stability, we use the fact that T (•) is simply monotone,

(λ -λ) T (x -y) ≥ 0, ∀x, y ∈ Ω, ∀λ ∈ T (x), ∀λ ∈ T (y). ( 59
)
in conjunction with the fact that there exists a ball of radius r ∈ (0, 4) included in T (0):

B r (0) = {z | z ≤ r} ⊂ T (0). (60)
For -λ ∈ T (z) and -λ 0 ∈ T (0), we get

(λ -λ 0 ) T z ≤ 0. (61)
Choosing z = 0 and -λ 0 = -rz/ z ∈ B r (0)) ⊂ T (0), we have

λ T z ≤ r z T z z = r z . (62)

Inria

Since -λ ∈ T (z), we get for z = 0,

V (z) + r |V (z)| ≤ 0. (63)
Assume for any t 1 > t 0 that V (z(t 1 )) > 0. Then we can divide (63) by |V (z)| and perform a time integration between t 1 and t 0

|V (z(t 1 ))| -|V (z(t 0 ))| ≤ -r(t 1 -t 0 ). (64)
From ( 64), we infer that lim t1→+∞ V (z(t 1 )) = -∞. By contradiction, we can conclude that there exists a finite time t 1 such that V (z(t 1 )) = 0. In other words, the equilibrium z = 0 is reached in finite time.

In the following section, we are interested in the discrete-time properties of the numerical scheme.

Properties of the numerical simulations

In Figure 4, several trajectories of the MCS [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF], which corresponds to the AP-extension of (4), have been simulated with the proposed algorithm. The simulations are been carried out over the time-interval [0, 3] with a time-step equal to h = 10 -2 . The integration parameter τ is chosen as τ = 1/2 (mid-point scheme).

Note that the equilibrium points are perfectly reproduced by the numerical simulation. Numerical simulations also reveals two sliding surfaces

{x 1 = θ 2 1 , θ 1 2 < x 2 } and {x 2 = θ 2 2 , θ 1 1 < x 1 }.
We can note that the simulated trajectories perfectly slide on these surfaces without any numerical chattering. As defined in [START_REF] Acary | Chattering-free digital sliding-mode control with state observer and disturbance rejection[END_REF], numerical chattering corresponds to the oscillations (limit cycles) which are solely due to the time-discretization of [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]. Avoiding numerical chattering is mainly due to the implicit character of the scheme for the inclusion part. Semi-explicit schemes of the form   

x k+1 = x k -h diag(γ)x k+τ + h g(σ k ) y k = Cx k -θ y k ∈ N [0,1] p (σ k ) (65)
yield numerical chattering. Although the convergence of such schemes is proven [START_REF] Dontchev | Difference methods for differential inclusions: a survey[END_REF], this drawback prevents qualitative analysis of solutions for a finite time-step size h > 0. Furthermore, the fact that the matrix -C B is definite positive in this area make it possible to apply Lemma 3 from [START_REF] Acary | Implicit Euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], which ensures that numerical chattering will not occur. In Figure 5, we show the simulation of the AP-extension of the PWL systems with the explicit Euler scheme (65) for τ = 0 (Figure 5(a)-(c)) and with the proposed implicit scheme (Figure 5(d)-(f)). With a 10 times larger time-step, the implicit scheme reaches the sliding surface and then the equilibrium without chattering. More precisely, there exists k 0 , such that x k = θ and (x k+1 -x k )/h = 0 for all k > k 0 .

In the reduced domain Ω, the one-step problem can be reduced to

x k+1 = x k -h diag(γ)x k+τ + h B σk+1 + hκ Cx k+1 -θ, ∈ N [0,1] 2 (σ k+1 ). ( 66
)
The Jacobian matrix is then given by The different trajectories illustrate that the system has three equilibrium points (of which are two stable and one unstable).

C B = -κ 1 0 0 -κ 2 . ( 67

Inria

Since κ i > 0, the matrix -C B is definite positive in Ω and we have

-x T C Bx ≥ κ x 2 (68)
with κ = min(κ 1 , κ 2 ). This means that the VI

θ -Cdiag(1/(1 + hτ γ)) (I n -h(1 -τ ) diag(γ))x k + h B σ + hκ) ∈ N [0,1] 2 (σ) (69)
is strongly monotone and has a unique solution (see [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF]Theorem 2.3.3]), given by To conclude this first example, the MCS formulation with an implicit time-discretization allows us to reproduce the main features of the dynamics of the PWL models of the gene reg-ulatory network: equilibrium points, sliding surfaces, finite-time stability. In the next section, some more realistic examples are considered.

(σ 2 1 ) k+1 = proj [0,1] 1 + hτ γ 1 hκ 1 θ 2 1 - 1 1 + hγ 1 (1 -h(1 -τ )γ 1 x 1,k+1 + hκ 1 ) (σ 2 2 ) k+1 = proj [0,1] 1 + hτ γ 2 hκ 2 θ 2 2 - 1 1 + hγ 2 (1 -h(1 -τ )γ 2 x 2,k+1 + hκ 2 ) . (70) 5 6 7 8 9 10 5 6 7 8 9 10 x2 x1 θ 2 2 θ 2 1 (a) 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 8.2 7.9 8 8.1 8.2 8.3 8.4 8.5 x2 x1 θ 2 2 θ 2 1 (b) 7.95 7

Applications in synthetic biology

In this section we apply the simulation methods developed above to the analysis of actual biological networks. We consider three examples from the field of synthetic biology [START_REF] Khalil | Synthetic biology: applications come of age[END_REF], which all concern relatively small networks that have been designed from pre-existing, well-characterized molecular elements and implemented in the cell to perform a particular function. The first two networks have been shown able to exhibit oscillations ( [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF][START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF], see [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF] for a review on synthetic oscillators). The third network is capable of a variety of dynamic behaviors and has been proposed as a benchmark for the reverse engineering of the network structure from time-course data [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF]. In the three cases we transformed the network structure into PWL equations, taking inspiration from existing ODE models, and we compared the dynamic features of the networks revealed by the simulations with available experimental data.

For each of the networks described above, we determined the equilibrium points and numerically computed solutions exemplifying the dynamics of the PWL systems. The mathematical analysis of equilibrium points of PWL models (1) was based on analytical results in the literature, concerning both equilibrium points within regions separated by threshold planes (x ∈ Θ) and equilibrium points located on one or several threshold planes (x ∈ Θ) [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF][START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF][START_REF] Plahte | Analysis and generic properties of gene regulatory networks with graded response functions[END_REF]. While the former are always asymptotically stable in the classical sense [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF], analyzing the stability of the latter requires generalized definitions adapted to the non-uniqueness of solutions of PWL systems. In particular, we used the results developed in [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] for F-extensions of PWL models (Definition 5), which provide easy-to-verify criteria for checking the stability in a number of typical cases. Notice that in the examples below, Assumptions 2 and 4 are always satisfied, which means that the F-extensions used for stability analysis and the AP-extensions used for simulation are equivalent (Proposition 11). As a consequence, the results obtained for the three example networks are valid for the different solution concepts reviewed in Section 3.

Repressilator

A repressilator is a network of several genes with a cyclic structure. Each gene codes for a transcription regulator that represses the next gene in the cycle. Elowitz and Leibler [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] have implemented a repressilator in the model bacterium E. coli, using the genes lacI, tetR, and cI (Figure 6(a)). The network can be externally controlled by adding the small inducer molecules IPTG and aTc to the growth medium, which inactivate LacI and TetR, respectively. Starting from their ODE model, we formulated a PWL model of the three-gene repressilator. The major change consists in lumping transcription and translation into a single step and replacing the sigmoidal regulation functions by step functions. The parameter values in [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] have been adapted accordingly (see Figure 6(b)).

The PWL model of the repressilator has a single equilibrium point, in accordance with the ODE model of Elowitz and Leibler. This equilibrium point is located at the intersection of the three thresholds, that is, at x = θ = [θ 1 , θ 2 , θ 3 ] T . The equilibrium point is thus also a discontinuity point of the PWL model. More precisely, the set of solutions starting from x = θ includes solutions that remain at the equilibrium point as well as solutions that leave this point at any time t ≥ 0. For the chosen parameter values, the equilibrium point is therefore unstable [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. The numerical solver finds an infinite number of solutions starting from x = θ. At each time-step, starting from the equilibrium point x = θ, an enumerative MCP solver is able to list

Inria LacI lacI TetR tetR cI CI (a)    ẋ1 = -γ1 x1 + κ 1 1 + κ 2 1 s -(x3, θ3) ẋ2 = -γ2 x2 + κ 1 2 + κ 2 2 s -(x1, θ1) ẋ3 = -γ3 x3 + κ 1 3 + κ 2 3 s -(x2, θ2) (71) (b)
Figure 6: (a) Repressilator consisting of three genes [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. The external control by the inducer molecules IPTG and aTc is not shown in the figure . (b) PWL model corresponding to the network in (a). The variables x 1 , x 2 , and x 3 represent the concentrations of the proteins LacI, TetR, and CI, respectively. The following parameter values have been used in the simulations:

γ 1 = γ 2 = γ 3 = 0.2, κ 1 1 = κ 1 2 = κ 1 3 = 4.8 10 -4 , κ 2 1 = κ 2 2 = κ 2 3 = 4.8 10 -1
, and θ 1 = θ 2 = θ 3 = 1. The time and concentration variables as well as the parameters in the model have been rescaled so as to make them dimensionless [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF].

multiple possible solutions which corresponds to staying at the equilibrium or leaving it. This qualitative feature of the algorithm will be illustrated in detail in Section 6.2. Figure 7 shows an example trajectory of the model, starting from the initial state x = 0 = [0, 0, 0] T . As can be seen, the trajectory approaches a limit cycle, again in accordance with the ODE model of Elowitz and Leibler. In-vivo measurements of the expression dynamics of a fluorescent reporter gene, driven by the same promoter as the gene cI, indicate that 40% of the analyzed cells indeed exhibit oscillations. The fact that not all cells oscillate, and that the periods and amplitudes vary between oscillating cells, can be attributed to stochastic effects that are not accounted for in deterministic ODE or PWL models [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF]. The occurrence of limit cycles in repressilator models has been an active subject of research for ODE models (see [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF][START_REF] Tyson | The dynamics of feedback control circuits in biochemical pathways[END_REF] for reviews) and PWL models [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Edwards | Control design for sustained oscillation in a two gene regulatory network[END_REF][START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: The case of a negative feedback loop[END_REF][START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF].

Synthetic oscillator with positive feedback

A second example of a synthetic oscillator network, developed by Atkinson et al. [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF], is shown in Figure 8(a). It consists of a negative feedback involving two genes, lacI and glnG: LacI inhibits the expression of glnG, whereas the product of the latter gene, the (phosphorylated) form of NRI (NRIp), activates lacI. The authors used a mutant strain which made NRI phosphorylation independent from the growth conditions, notably the cellular nitrogen status. NRIp also activates the expression of its own gene, thus giving rise to an additional positive feedback loop. While a negative feedback loop is necessary for generating oscillations, a positive feedback loop is believed to favor the robustness of the oscillations [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF]. Like in the case of the repressilator, the activity of LacI can be externally controlled by varying the concentration of IPTG in the growth medium.

We constructed a PWL model of this network, given by (72) in Figure 8(b). The model has been obtained by first formulating an ODE model using information from the paper, in the same way as for the repressilator. In particular, the model includes the design features that LacI and NRIp have antagonistic effects on glnG and that the binding of LacI to the promoter region shuts off the expression of the gene. Moreover, the affinity of NRIp for the promoter region upstream of glnG is higher than for the promoter region upstream of lacI. This gives rise to a lower value of Figure 7: Simulation of repressilator network in Figure 6: (a) a solution starting from the initial state x 0 = 0 and (b) the corresponding phase-space trajectory. Notice that time is a dimensionless variable in the model [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. Physical time can be reconstructed by multiplying the rescaled time variable in the model by 2/ ln 2. The dashed lines represent the threshold concentrations of the variables.

Inria the threshold associated with NRIp-activation of glnG than the threshold associated with NRIpactivation of lacI (θ 1 2 < θ 2 2 ). In a second step, this ODE model has been simplified by lumping the transcription and translation steps, and by replacing the sigmoidal regulation functions by step functions. The parameter values were copied from the model in [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF], and where necessary completed with reasonable estimates.

LacI lacI NRI glnG NRIp (a) ẋ1 = -γ1 x1 + κ1 s + (x2, θ 2 2 ) ẋ2 = -γ2 x2 + κ2 s -(x1, θ1) s + (x2, θ 1 2 ) (72) (b)
Figure 8: (a) Oscillator with positive feedback consisting of two genes [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF]. The external control by the inducer molecule IPTG is not shown in the figure . (b) PWL model corresponding to the network in (a). The variables x 1 and x 2 represent the concentrations of the proteins LacI and GlnG, respectively. The following parameter values have been used in the simulations: γ 1 = γ 2 = 0.032, κ 1 = 0.08, κ 2 = 0.16, and θ 1 = 1, θ 1 2 = 1, and θ 2 2 = 4. The time and concentration variables as well as the parameters in the model have been rescaled so as to make them dimensionless, in the same way as in [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. The dashed lines represent the threshold concentrations of the variables.

The PWL model (72) has three equilibrium points, with the coordinates (0, 0) T , (0, θ 1 2 ) T , and (θ 1 , θ 2 2 ) T , respectively. The local dynamics around these equilibrium points are shown in Figure 9(a). The first equilibrium point lies away from the threshold planes and is asymptotically stable [START_REF] Glass | The logical analysis of continuous non-linear biochemical control networks[END_REF]. The second and the third equilibrium points, on the other hand, are located on a threshold plane and on the intersection of two threshold planes, respectively. More precisely, the equilibrium point x = (0, θ 1 2 ) T lies on a repulsive discontinuity segment, defined by 0 ≤ x 1 < θ 1 and x 2 = θ 1 2 . As a consequence, it is unstable [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. Figure 9(a) shows that solutions starting in the neighborhood of x = (θ 1 , θ 2 2 ) T converge towards this point. This illustrates that, for the chosen parameter values, this equilibrium point is asymptotically stable.

The existence of a repulsive discontinuity segment in Figure 9 entails non-uniqueness of solutions. In fact, any solution with an initial state x 0 in the segment 0 < x 1 < θ 1 and x 2 = θ 1 2 has an infinite number of solutions. The set of solutions includes one that, while sliding on the segment, asymptotically converges towards the equilibrium point x = (0, θ 1 2 ) T . The other solutions all slide for some time ℵ ≥ 0 on the segment and then escape to either the region above or below. That is, they converge towards one of the two asymptotically stable equilibrium points of the system, x = (0, 0

) T or x = (θ 1 , θ 2
2 ) T ). The discontinuity segment is the PWL analog of a separatrix in classical ODE models. It occurs as a consequence of the positive autoregulatory interaction in the circuit of Figure 8, corresponding to the step function s + (x 2 , θ 1 2 ) in the righthand side of (72).

In Figure 9(a), phase-space trajectories corresponding to several solutions starting from the initial state x 0 = (0.78, 1) T located on the repulsive discontinuity segment are depicted. In practice, we first performed a simulation with an enumerative solver that reveals that in the neighborhood of the repulsive segment the MLCP (49) has multiple solutions. The modes, as defined in Section 4, are recorded. Once this first simulation has been performed, the user can choose a preferred solution (encoded by one triplet of index sets) on a user-supplied part of the trajectory, by asking the numerical solver to check if a solution with the given triplet exists. For the trajectories shown in Figure 9(a), three possible solutions have been detected in the neighborhood of the repulsive segment. A new simulation has been performed asking to prefer the triplet corresponding to the solution that slides on the segment. For the other trajectories, we changed the preferred solution when the trajectory reaches x 1 = 0.55 and x 1 = 0.4 to the triplet of solutions that yields a trajectory escaping below and above, respectively.

Atkinson et al. [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF] measured the time-varying dynamics of the circuit on the population level using the chromosomal lacZ gene. This gene is under the control of LacI and encodes β-galactosidase, whose levels can be assayed. The authors observed damped oscillations of βgalactosidase, which corresponds to the dynamics of the PWL model when the initial LacI concentration lies above its first threshold (x 0 1 > θ 1 1 , see Figure 9). In order to verify that damped oscillations on the population level are not due to the desynchronization of individual cells that each exhibit stable oscillations, Atkinson et al. [START_REF] Atkinson | Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli[END_REF] also measured gene expression in single cells, using a fluorescent reporter gene fused to a LacI-repressible promoter. The results confirm the occurrence of damped oscillations in the oscillator with positive feedback, in agreement with the population-level assay.

IRMA: a synthetic benchmark network

IRMA is a synthetic network constructed in yeast and proposed as a benchmark for modeling and identification approaches [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF]. The network consists of five well-characterized genes that have been chosen so as to include different kinds of interactions, notably transcription regulation and protein-protein interactions. The structure of the IRMA network is shown in Figure 10(a). The expression of the CBF1 gene is positively regulated by Swi5 and negatively regulated by Ash1. CBF1 encodes the transcription factor Cbf1 that activates expression of the GAL4 gene. The expression of SWI5 is activated by Gal4, but only in the absence of Gal80 or in the presence of galactose. Gal80 binds to the Gal4 activation domain, but galactose releases this inhibition of transcription. The product of SWI5 activates the expression of three other genes in the network: ASH1, CBF1, and GAL80. The network contains one positive and two negative feedback loops. Consequently, for suitable parameter values IRMA might function as a synthetic oscillator [START_REF] Marucci | How to turn a genetic circuit into a synthetic tunable oscillator, or a bistable switch[END_REF].

Figure 10(b) shows a PWL model of IRMA that was previously developed [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF]. The regulatory functions in the model encode the interaction structure described above. Most parameters values in the PWL model were taken from the original IRMA publication. In a few cases, however, they were slightly adapted to satisfy inequality constraints between parameters that were inferred from the experimental data [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF].

Using the model (73), we simulated the behavior of the network in response to two different perturbations: shifting cells from glucose to galactose medium (switch-on experiments), and from galactose to glucose medium (switch-off experiments). The terms switch-on (switch-off) refer to the activation (inhibition) of SWI5 expression during growth on galactose (glucose). For these two perturbations, the temporal evolution of the expression of all the genes in the network was monitored by qRT-PCR with good time resolution in [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF].

In switch-off conditions, when s + (u, θ u ) = 1, the PWL system has a single asymptotically stable equilibrium point, not located on a threshold plane. At this equilibrium point all genes are expressed at their basal level x = (κ 1 1 /γ 1 , κ 1 2 /γ 2 , κ 1 3 /γ 3 , 0, κ 1 5 /γ 5 ) T . The measured time-series expression profiles of the five network genes in a switch-off experiments indeed show that the system evolves towards a basal state [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF].

In switch-on conditions, when s + (u, θ u ) = 0, the PWL system has three equilibrium points.

2 ) T = (1, 4) T . Notice that time is a dimensionless variable in the model [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. Physical time can be reconstructed by multiplying the rescaled time variable in the model by 2/ ln 2.

Discussion

We have presented a numerical method for the simulation of gene regulatory networks described by PWL models. In order to account for discontinuities in the right-hand side of the differential equations, PWL models have to be extended to differential inclusions which cannot be handled by classical ODE solvers. We show that by reformulating one such extension, the AP-extension, as a mixed complementarity system (MCS), numerical methods developed for the time-integration of nonsmooth dynamical systems are directly applicable. An extensive amount of work on the mathematical and numerical analysis of complementarity systems has been done, which provides a solid theoretical basis for the proposed method. Moreover, powerful and widely-used numerical tools for complementarity systems are available, like the Siconos platform used here.

We have shown that under assumptions that are not restrictive in practice, the AP-extensions are equivalent with other proposed extensions of PWL models, notably classical F-extensions (Proposition 11). The existence of solutions, which is guaranteed for F-extensions, therefore carries over to AP-extensions. Under the same assumptions, the AP-extensions are also equivalent with hyperrectangular overapproximations of F-extensions that have been used in qualitative simulations (Corollary 12). This has the advantage that the numerical analysis of PWL models can be preceded by the symbolic computation of equilibrium points and (in some cases) the analytic determination of their stability [START_REF] Jong | Search for steady states of piecewise-linear differential equation models of genetic regulatory networks[END_REF], as was done for the examples discussed in Section 6. More generally, the equivalency of the above-mentioned differential inclusions shows that the numerical method developed here is applicable for different solution concepts of PWL models proposed in the literature.

The practical usefulness of the approach was illustrated by means of the numerical analysis of three actual biological networks in synthetic biology. We developed PWL models for the circuits, either from scratch or by adapting existing ODE models, and parametrized the models with the help of literature data. For each of the models, we determined the equilibrium points and their stability, and we numerically computed solutions of the characteristic time-behavior of the networks. Although PWL models are coarse-grained approximations of the actual biological processes taking place in the cell, and the dynamics of PWL models may differ from those of ODE models [START_REF] Polynikis | Comparing different ODE modelling approaches for gene regulatory networks[END_REF], many of the observed dynamical properties of the synthetic networks could be reproduced by the simulations, such as the occurrence of stable oscillations in the repressilator (Section 6.1) and damped oscillations in the oscillator with positive feedback (Section 6.2).

The formulation of PWL models in the complementarity systems framework provides access to powerful numerical tools for the analysis of dynamic properties of naturally-occurring networks [START_REF] Abou-Jaoude | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF][START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF] or for the design of control schemes for synthetic networks [START_REF] Edwards | Control design for sustained oscillation in a two gene regulatory network[END_REF]. Moreover, it also opens up the possibility to exploit results obtained for complementarity systems in the context of the analysis of PWL models. For example, as a first step in this direction, we showed in Section 5.1 that the strong monotonicity of the differential inclusion obtained after reformulating the PWL model as a MCS can be used to assure that an equilibrium point is (finite-time) stable.

A Alternative AP-extension of PWL models

As discussed in Section 3.1, the definition of AP-extensions could be relaxed by substituting different values σ k j ∈ S + (x j , θ k j ) for the different occurrences of s + (x j , θ k j ), and by decoupling the values for positive and negative step functions. This gives rise to an alternative AP-extension of the PWL model [START_REF] Abou-Jaoude | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF].

Definition 18 (Alternative AP-extension of PWL models). The alternative AP-extension of a

Inria

Figure 1 :

 1 Figure 1: (a) Example of a gene regulatory network of two genes, each coding for a regulatory protein. (b) PWL model corresponding to the network in (a).

Proposition 3 .

 3 All PWL models (1) with regulation functions b l i (•) in DNF satisfy Assumption 2.

Figure 2 (Figure 2 :

 22 Figure 2: The F-and AP-extensions of the PWL model (18) are different at x = (θ 1 , θ 2 ) T , the intersection point of the thresholds. (a) F(x), defined by (19), is the smallest closed convex set containing the vector fields of the four neighboring regions of x = (θ 1 , θ 2 ) T , evaluated at the intersection point. (b) G(x) is obtained by varying the parameters σ 1 , σ 2 over the interval [0, 1],as defined in[START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF]. It includes the vector fields in the neighboring regions of x = (θ 1 , θ 2 ) T , which are obtained when σ 1 , σ 2 take their extreme values 0 or 1. The plots are obtained for κ 1 > γ 1 θ 1 and κ 2 > γ 2 θ 2 .

Figure 3 :

 3 Figure 3: Right-hand side of the Filippov differential inclusion (21) for the example network of Figure 1, at x = (x 1 , x 2 ) on the segment {θ 1 1 < x 1 < θ 2 1 , x 2 = θ 2 2 }. In this example, F(x) and G(x) coincide.

Figure 4 :

 4 Figure 4: Different trajectories of the two-gene regulatory network in Figure 1 defined by (4). The different trajectories illustrate that the system has three equilibrium points (of which are two stable and one unstable).

Figure 5 :

 5 Figure 5: Illustration of the numerical chattering effect by zooming in on the state space on the system (4) as shown in Figure 4. The trajectories start from x(t 0 ) = [10, 5] T . (a)-(c). Different zooms on the trajectory obtained with the explicit Euler scheme (h = 10 -3 ). (d)-(f). Different zooms on the trajectory obtained with the proposed implicit Euler scheme (h = 10 -2 ).

Figure 9 :

 9 Figure 9: Simulation of synthetic oscillator network with positive feedback in Figure 8. (a)Trajectories corresponding to several solutions starting from the initial state x 0 = (0.78, 1) T located on the repulsive discontinuity segment separating the basins of attractions of the two asymptotically stable equilibrium points. (b) Simulation trace of the damped oscillation that corresponds to the trajectory converging towards x = (θ 1 , θ2 2 ) T = (1, 4) T . Notice that time is a dimensionless variable in the model[START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. Physical time can be reconstructed by multiplying the rescaled time variable in the model by 2/ ln 2.

RR n°8207

The regulatory logic can be inferred by writing the equations in DNF. In the case of b 1 (x) = s + (x 1 , θ 1 ) + s + (x

, θ 2 ) -2s + (x 1 , θ 1 ) s + (x 2 , θ 2 ) this yields s + (x 1 , θ 1 ) s -(x 2 , θ 2 ) + s -(x 1 , θ 1 ) s + (x 2 , θ 2 ), which corresponds to a Boolean XOR function. The function b 2 (x) = 1 -s + (x 1 , θ 1 ) s + (x 2 , θ 2 ) is the algebraic equivalent of a Boolean NAND, that is, s -(x 1 , θ 1 ) s -(x 2 , θ 2 ) + s + (x 1 , θ 1 ) s -(x 2 , θ 2 ) + s -(x 1 , θ 1 ) s + (x 2 , θ 2 ). RR n°8207

http://siconos.gforge.inria.fr Inria

Figure 10: (a) IRMA network consisting of five genes [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF]. (b) PWL model corresponding to the network in (a). The variables x 1 , x 2 , x 3 , x 4 , and x 5 represent the concentrations of the proteins Gal4, Swi5, Ash1, Cbf1, and Gal80, respectively. Like in [START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering[END_REF], these variables have arbitrary units (a.u.). In addition, the variable u represents the concentrations of external galactose, an input signal. The following parameter values have been used in the simulations: γ 1 = 0.05, γ 2 = 0.04, γ 3 = 0.05, γ 4 = 0.02, γ 5 = 0.6, κ 1 1 = 1.1 10 -4 , κ 2 1 = 9.0 10 -4 , κ 1 2 = 3.0 10 -4 , κ 2 2 = 0.15, κ 1 3 = 6.0 10 -4 , κ 2 3 = 0.018, κ 1 4 = 5.0 10 -4 , κ 2 4 = 0.03, κ 1 5 = 7.5 10 -3 , κ 2 5 = 0.015, θ 1 = 0.01, θ 1 2 = 0.01, θ 2 2 = 0.06, θ 3 2 = 0.08, θ 3 = 0.035, θ 4 = 0.04, and θ 5 = 0.01. The units of the degradation parameters γ, κ, and θ are min -1 , a.u., and a.u. min -1 , respectively. We do not assign numerical values to u and its threshold θ u , but rather consider separately the cases s + (u, θ u ) = 0 (switch-on) and s + (u, θ u ) = 1 (switch-off). The first coincides with the equilibrium point in switch-off conditions. The second is an unstable equilibrium point located on the intersection of three different threshold planes

The third equilibrium point, with the coordinates (θ 1 , θ 2 2 , θ 3 , θ 4 , κ 1 5 /γ 5 ) T is also unstable: solutions starting in its neighborhood converge towards a limit cycle.

Figure 11 shows the simulated response of the network in a switch-on experiment, with initial conditions that are obtained by slightly perturbing the switch-off equilibrium point. As can be seen, the solutions do not contain sliding modes, in which one or more variables remain at their threshold value during a finite-time interval. The agreement of the predicted oscillations with the experimental data is unclear. Individual time-series, which were used for the parameterization of the PWL model [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF], seem to exhibit oscillatory patterns. These patterns largely disappear when the data are averaged over several independent time-series (where damped oscillations are seen for Swi5 and Cbf1 only). The large error bars in the averaged experiments may hide more subtle dynamic patterns, however, no data on the single-cell level are available.

PWL model (1) is defined by the following differential inclusion

. . .

Ĝn (x)

where S = (S + (x 1 , θ 1 1 ), . . . , S + (x 1 , θ p1 1 ), . . . , S + (x n , θ 1 n ), . . . , S + (x n , θ pn n )) T , and bl i (•) is obtained from b l i (•), i ∈ {1, . . . , n}, by replacing every occurrence of s + (x j , θ k j ) and s -(x j , θ k j ) by S + (x j , θ k j ) and S -(x j , θ k j ), respectively. This definition seems similar to Definition 8 but not equivalent, as shown in the example below. Consider the following simple PWL model:

and let us assume that

The model (75) has three distinct regimes:

1.

, the dynamical equation is given by: ẋ1

The dynamical equation is given by: ẋ1

3. θ 1 < x 1 . The dynamical equation is given by:

Notice that x 1 = θ 1 is a discontinuity point of the system. In order to model the dynamics at this point, first consider the differential inclusion given by Definition 8:

with

Under the assumption (76), it therefore holds that ẋ1 < 0, as expected from the dynamics in the regions above and below x 1 = θ 1 (see Figure 12). The point of discontinuity is a crossing point in the terminology of Filippov solutions. A quick analysis reveals that, when the differential inclusion is defined as in Definition 8, the system has a single equilibrium point x 1 = κ 1 1 /γ 1 . Alternatively, following (74), the dynamics of the PWL model can be modeled by means of the differential inclusion

At

Note that (83) differs from (81), and that 0 ∈ Ĝ(x 1 ) since κ 1 1 + κ 2 1 > γ 1 θ 1 . This means that the differential inclusion (83) introduces an additional equilibrium point in θ 1 which is not a solution in the sense of Definition 8 (Figure 12). As explained in Sec. 3.1, we prefer Definition 8 above Definition 18.