Camilo La Rota
email: camilo.larota@ens-lyon.fr

Fabien Tarissan
email: tarissan@lix.polytechnique.fr

Leo Liberti
email: liberti@lix.polytechnique.fr

) -Book of Extended Abstracts INFERRING PARAMETERS IN GENETIC REGULATORY NETWORKS

INTRODUCTION

We formulate and solve the problem of determining the arc weights in a Genetic Regulatory Network (GRN) so that it presents observed states at given development stages. Complex biological regulatory networks determine the function and some observed emergent phenomena in most biological systems. In particular, molecular regulatory networks (proteins, mR-NAs, hormones, etc) underlie morphogenesis phenomena. It is therefore very important to understand the dynamics of these networks in order to infer the molecular mechanisms regulating plant and animal development. In the present work, we use experimental data to reconstruct the GRN acting in each cell. More precisely, we attempt to determine some possible scenarios of development together with candidate gene sub-networks that we assume can partially explain the expression patterns observed at each development stage. The estimated values of the kinetic parameters (arc weights) must be such that the stable states (fixed points with respect to the network evolution dynamics) of each sub-network reflect the spatial patterns of gene expression (molecule presence or absence) observed at each development stage. In particular, we apply this approach to the inference of the regulatory network controlling the early stages of carpel development in the model plant Arabidopsis thaliana and we show some results of this work. This network is still poorly understood, even if data concerning genes and morphogenes involved, single regulatory interactions and small circuits exist in the literature [START_REF] Balanzá | Patterning the female side of arabidopsis: the importance of hormones[END_REF].

MODELLING THE PROBLEM

Given a directed graph G = (V, A), a discrete set of time instants T (which we suppose to be an initial contiguous proper subset of N) and the following functions:

• a function α : A → {+1, -1} called the arc sign function; • a function ω : A → R + called the arc weight function; • a function x : V × T → {0, 1} called the gene activation function; • a function ι : V → {0, 1} called the initial configuration; • a function θ : V → R + called the threshold function, A gene regulatory network (GRN) is a 7-tuple (G, T, α, ω, x, ι, θ) such that: ∀v ∈ V x(v, 1) = ι(v) (1)
∀v ∈ V, t ∈ T {1} x(v, t) =    1 if u∈δ -(v) α(u, v)ω(u, v)x(u, t -1) ≥ θ(v) 0 otherwise, (2)
where 1)-( 2) together are called the evolution rules of the GRN. For any particular t ∈ T , x(•, t) : V → {0, 1} is called a configuration.

δ -(v) = {u ∈ V | (u, v) ∈ A} for all v ∈ V . Eqns. (
Since the evolution rules relate a configuration at time t with a configuration at time t -1, if x(•, t) = x(•, t -1) then x(•, t ′ ) = x(•, t) for all t ′ > t: such configurations are called fixed points of the GRN.

In this paper we deal with an inverse problem relating to the reconstruction of a stable subnetwork in GRNs from observed data. More precisely, we address the following.

Stable Subnetwork Reconstruction in GRNs (SSRGRN). Given a digraph G = (V, A), a time instant set T , an arc sign function α, an initial configuration C 0 : V → {0, 1}, a set U ⊆ V determining the nodes of the (induced) subnetwork of G to be reconstructed and the stable configurations S : U → {0, 1} sought for the subnetwork, find an arc weight function ω and a threshold function θ with the property that for all ι ∈ C 0 there exists a gene activation function x such that (G, T, α, ω, x, ι, θ) is a GRN whose fixed points are at a minimum distance to observed data.

In other words, we attempt to estimate the arc weights and threshold functions of a GRN from the knowledge of the digraph topology G, the induced subnetwork topology U and the arc sign function α in such a way that (a) the GRN evolution rules are consistent with respect to a certain set of initial configurations and (b) the fixed points in the subnetwork induced by the estimated values are as close as possible to the given ones. The methodology we shall follow is that of modelling the SSRGRN by means of a mathematical programming formulation:

min x f (x) subject to g(x) ≤ 0, (3)
where x ∈ R n are the decision variables and f : R n → R is the objective function to be minimized subject to a set of constraints g : R n → R m which may also include variable ranges or integrality constraints on the variables. As shown below, the SSRGRN is a Mixed-Integer Nonlinear Program (MINLP) which can be reformulated to a Mixed-Integer Linear Program (MILP). Our primary concern in solving the SSRGRN is thus modellistic rather than algorithmic. One of the foremost difficulties is that of employing a static modelling paradigm -such as mathematical programming -in order to describe a problem whose very definition depends on time. Another important difficulty resides in expressing the necessary and sufficient conditions for a configuration to be a fixed point as constraints suitable for use in a formulation like (3). In the following, we present a simplified version in which only one targeted configuration S is used.

ACKNOWLEDGMENTS

This work was supported by EU FP6 FET Project "Morphex". We are also deeply grateful to Monéger and Traas from the RDP-ENS-Lyon institute for having provided the data on which we tested our approach.

2. fixed point conditions:

REFORMULATION AND SOLUTION

The above problem is a nonconvex MINLP that can be reformulated exactly to a MILP using the techniques proposed in [START_REF] Liberti | Reformulation techniques in mathematical programming[END_REF]. We solved to optimality a few real-life instances from the GRN of Arabidopsis thaliana using AMPL [START_REF] Fourer | The AMPL Book[END_REF] to model the problem and CPLEX [START_REF] Ilog | ILOG CPLEX 10.1 User's Manual[END_REF] to solve it. The size of the GRNs involved were such that CPLEX obtained the optimal solution in a matter of seconds and/or minutes.


Nicolas Champagnat
email: nicolas.champagnat@inria.fr

Denis Villemonais
email: denis.villemonais@univ-lorraine.fr

Exponential convergence to quasi-stationary distribution for absorbed one-dimensional diffusions with killing

Keywords: diffusions, one-dimensional diffusions with killing, absorbed process, quasi-stationary distribution, Q-process, uniform exponential mixing property, one dimensional processes with jumps. 2010 Mathematics Subject Classification. Primary: 60J60, 60J70, 37A25, 60B10, 60F99. Secondary: 60G44, 60J75

This article studies the quasi-stationary behaviour of absorbed onedimensional diffusion processes with killing on [0, ∞). We obtain criteria for the exponential convergence to a unique quasi-stationary distribution in total variation, uniformly with respect to the initial distribution. Our approach is based on probabilistic and coupling methods, contrary to the classical approach based on spectral theory results. Our general criteria apply in the case where ∞ is entrance and 0 either regular or exit, and are proved to be satisfied under several explicit assumptions expressed only in terms of the speed and killing measures. We also obtain exponential ergodicity results on the Q-process. We provide several examples and extensions, including diffusions with singular speed and killing measures, general models of population dynamics, drifted Brownian motions and some one-dimensional processes with jumps.

Introduction

This article studies the quasi-stationary behaviour of general one-dimensional diffusion processes with killing in an interval E of R, absorbed at its finite boundaries. When the process is killed or absorbed, it is sent to some cemetary point ∂. This covers the case of solutions to one-dimensional stochastic differential equations (SDE) with space-dependent killing rate, but also of diffusions with singular speed and killing measures.

We recall that a quasi-stationary distribution for a continuous-time Markov process (X t , t ≥ 0) on the state space E ∪ {∂}, is a probability measure α on E such that

P α (X t ∈ • | t < τ ∂ ) = α(•), ∀t ≥ 0,
where P α denotes the distribution of the process X given that X 0 has distribution α, and τ ∂ := inf{t ≥ 0 : X t = ∂}.

We refer to [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF][START_REF] Van Doorn | Quasi-stationary distributions for discrete-state models[END_REF][START_REF] Pollett | Quasi-stationary distributions: A bibliography[END_REF] for general introductions to the topic. Our goal is to give conditions ensuring the existence of a unique quasilimiting distribution α on E, i.e. a probability measure α such that for all probability measures µ on E and all A ⊂ E measurable,

lim t→+∞ P µ (X t ∈ A | t < τ ∂ ) = α(A), (1.1)
where, in addition, the convergence is exponential and uniform with respect to µ and A. In particular, α is the unique quasi-stationary distribution. This topic has been extensively studied for one-dimensional diffusions without killing in [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF][START_REF] Littin | Uniqueness of quasistationary distributions and discrete spectra when ∞ is an entrance boundary and 0 is singular[END_REF][START_REF] Miura | Ultracontractivity for Markov semigroups and quasistationary distributions[END_REF][START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF], where nearly optimal criteria are obtained. The case with killing is more complex and the existing results cover less general situations [START_REF] Steinsaltz | Quasistationary distributions for one-dimensional diffusions with killing[END_REF][START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]. In particular, these references are restricted to the study of solutions to SDEs with continuous absorption rate up to the boundary of E, and let the questions of uniqueness of the quasi-stationary distribution and of convergence in (1.1) open.

The present paper is focused on the case of a diffusion on E = [0, +∞) absorbed at 0 with scale function s, speed measure m and killing measure k, assuming that killing corresponds to an immediate jump to ∂ = 0. We consider the situation where ∞ is an entrance boundary and 0 is either exit or regular. Our results easily extend to cases of bounded intervals with reachable boundaries.

We give two criteria, each of them involving one condition concerning the diffusion (without killing) with scale function s and speed measure m, and another condition on the killing time. The condition on the diffusion without killing comes from [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF] and covers nearly all one-dimensional diffusions on [0, ∞) such that ∞ is an entrance boundary and a.s. absorbed in finite time at 0. The conditions on the killing time only concern the behavior of the diffusion and of the killing measure in the neighborhood of 0, as soon as ∞ is an entrance boundary. In order to apply these results to practical situations, we provide several explicit criteria ensuring all these conditions and a series of examples that enter our setting.

This contribution improves known results in several directions. First, it covers situations of diffusion processes which are not solutions to SDEs and cases of irregular and unbounded killing rates. Secondly, aside from proving the existence of a quasi-stationary distribution, we also obtain exponential convergence of conditional distributions. In particular, all the initial distributions of X belong to the domain of attraction of the quasi-stationary distribution. This is of great importance in applications, since in practice, one usually has no precise estimate of the initial distribution. In addition, our estimates of convergence are exponential and uniform in total variation norm, and hence provide a uniform bound for the time needed to observe stabilisation of the conditional distribution of the process, regardless of the initial distribution (see also the discussion in [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF]Ex. 2]).

The main methodological novelty of our proofs relies on its purely probabilistic approach. We do not use any spectral theoretical result (typically, Sturm-Liouville theory), which are the key tool of all the previously cited works on diffusions, except [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]. Instead, we use criteria for general Markov processes proved in [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF], based on coupling and Dobrushin coefficient techniques. These criteria also imply exponential ergodicity for the Q-process, defined as the process X conditioned to never be absorbed. The generality and flexibility of this approach also allow to cover, without substantial modification of the arguments, many situations where the spectral theory has received much less attention, such as one-dimensional diffusions with jumps and killing.

For our study, we need to give a probabilistic formulation of the property that ∞ is entrance. In the case without killing, this is known as the classical property of coming down from infinity (see e.g. [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF][START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]). For a diffusion with killing, we show that a process with entrance boundary at ∞ comes down from infinity before killing in the sense that the diffusion started from +∞ hits 0 in finite time before killing with positive probability.

The paper is organized as follows. In Section 2, we precisely define the absorbed diffusion processes with killing under study and give their construction as time-changed Brownian motions. Although quite natural, this alternative construction of killed diffusions is not given in classical references. Section 3.1 contains the statement of our main results on the exponential convergence of conditional distributions, the asymptotic behavior of the probability of survival and the existence and ergodicity of the Q-process. In Subsection 3.2, we give explicit criteria ensuring the conditions of Section 3.1. We focus in Sections 2 and 3 on diffusions on natural scale. The extension to general diffusions is given in Section 4.1, together with a series of examples of diffusions which are not solutions to SDEs, of diffusions with unbounded or irregular killing rates and of one-dimensional processes with killing and jumps. Section 4.2 is devoted to the study of models of population dynamics of the form

dY t = Y t dB t + Y t h(Y t )dt
and Section 4.3 to drifted Brownian motions, which are the basic models of [START_REF] Steinsaltz | Quasistationary distributions for one-dimensional diffusions with killing[END_REF][START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]. Section 5 concerns the property of coming down from infinity before killing. Finally, Section 6 gives the proofs of our main results.

Absorbed diffusion processes with killing

Our goal is to construct diffusion processes with killing on [0, +∞), absorbed at ∂ = 0. The typical situation corresponds to stochastic population dynamics of continuous densities with possible continuous or sudden extinction.

We recall that a stochastic process (X t , t ≥ 0) on [0, +∞) is called a diffusion (without killing) if it has a.s. continuous paths, satisfies the strong Markov property and is regular. By regular, we mean that for all x ∈ (0, ∞) and y ∈ [0, ∞), P x (T y < ∞) > 0, where T y is the first hitting time of y by the process X. Given such a process, there exists a continuous and strictly increasing function s on [0, ∞), called the scale function, such that (s(X t∧T 0 ), t ≥ 0) is a local martingale [START_REF] Freedman | Brownian motion and diffusion[END_REF]. The stochastic process (s(X t ), t ≥ 0) is itself a diffusion process with identity scale function. Replacing (X t , t ≥ 0) by (s(X t ), t ≥ 0), we can assume without loss of generality that s(x) = x.

To such a process X on natural scale, one can associate a unique locally finite positive measure m(dx) on (0, ∞), called the speed measure of X, which gives positive mass to any open subset of (0, +∞) and such that X t = B σt for all t ≥ 0 for some standard Brownian motion B, where

σ t = inf {s > 0 : A s > t} , with A s = ∞ 0 L x s m(dx) (2.1)
and L x is the local time of B at level x. Conversely, given any positive locally finite measure m on (0, ∞) giving positive mass to any open subset of (0, ∞), any such time change of a Brownian motion defines a regular diffusion on [0, ∞) [9, Thm. 23.9]. Note that, since σ t is continuous and since X t = B σt for all t ≥ 0 a.s., we have

σ T 0 = T B 0 , or, equivalently, T 0 = A T B 0 , (2.2)
where T B x is the first hitting time of x ∈ R by the process B. In this work, we study diffusion processes on [0, +∞) with killing as defined in [START_REF] Itô | Diffusion processes and their sample paths[END_REF]. As above, we can assume without loss of generality that the diffusion is on natural scale. We show below that such a process can be obtained from an explicit pathwise construction from a given Brownian motion (B t , t ≥ 0) and an independent exponential random variable E of parameter 1. Although this construction is quite natural, this is not done in the classical references [START_REF] Itô | Diffusion processes and their sample paths[END_REF][START_REF] Freedman | Brownian motion and diffusion[END_REF][START_REF] Méléard | Application du calcul stochastique à l'etude de processus de Markov réguliers sur [0, 1[END_REF][START_REF] Kallenberg | Foundations of modern probability. Probability and its Applications[END_REF].

Let k and m be two positive locally finite measures on (0, +∞), such that m gives positive mass to any open subset of (0, +∞). The measures m and k will be referred to as the speed measure and the killing measure, respectively, of the diffusion with killing.

We first define the diffusion process on natural scale without killing as above by

X t := B σt , ∀t ≥ 0,
where σ t is defined by (2.1). Next, we define for all t ≥ 0

κ t := ∞ 0 L y σt k(dy), ∀t ≥ 0. (2.
3)

It follows from the definition of the local time (L X,x t , t ≥ 0, x ≥ 0) of the process X of [6, p. 160] that

κ t = ∞ 0 L X,y t k(dy).
Let E be an exponential r.v. of parameter 1 independent of B. Then we define τ κ := inf{t ≥ 0 : κ t ≥ E} and T 0 := inf{t ≥ 0 : X t = 0}.

(2.4)

Our goal is to construct a process X that can hit 0 either continuously following the path of X, or discontinuously when it is killed at time τ κ . This leads to the following definition of the time of discontinuous absorption (or killing time) τ d ∂ and the time of continuous absorption τ c ∂ :

τ d ∂ := τ κ if τ κ < T 0 , +∞ if T 0 ≤ τ κ and τ c ∂ := +∞ if τ κ < T 0 , T 0 if T 0 ≤ τ κ . (2.5)
The absorption time τ ∂ is then defined as

τ ∂ := τ d ∂ ∧ τ c ∂ .
Then, the absorbed diffusion process with killing X is defined as

X t := X t if t < τ ∂ , ∂ = 0 otherwise.
(2.6)

Proposition 2.1. The process (X t , t ≥ 0) defined in (2.6) is a diffusion process with killing on natural scale, with speed measure m and killing measure k, as defined in [START_REF] Itô | Diffusion processes and their sample paths[END_REF].

Proof. Note that

κ t = ∞ 0 L X,y t k(dy), ∀t < τ κ ,
since X and X coincide before τ κ . Since E is independent of B,

P(τ κ > t | X t , t ≥ 0) = P(τ κ > t | B s , s ≥ 0) = exp - ∞ 0 L X,y t k(dy) .
(2.7) This is exactly the definition of the distribution of the killing time of the diffusion process with killing of [8, p. 179].

In particular, in the case where k is absolutely continuous with respect to m, the killing rate of the diffusion at position x (in the sense e.g. of [START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]) is given by dk dm (x). This can be deduced from (2.7) and from the occupation time formula of [START_REF] Freedman | Brownian motion and diffusion[END_REF].

In addition, the property of regularity of the diffusion with killing can be easily checked from this definition, as shown in the next proposition. Proposition 2.2. For all x, y > 0,

P x (T y < ∞) > 0,
and for all x > 0 and t > 0,

P x (τ ∂ > t) > 0, ∀x > 0, ∀t ≥ 0.
(2.8)

Proof. Fix x = y. On the event { T y < ∞}, we have σ Ty = T B y < T B 0 . In particular, the infimum of B s over [0, T B y ] is positive, and its maximum is finite. Hence the function y → L y T B y is a.s. continuous with compact support in (0, +∞), hence a.s. bounded, on the event { T y < ∞}. Since k is locally finite on (0, +∞) and P x ( T y < ∞) > 0, this implies that

P x (T y < ∞) = P x ( T y < ∞ and κ Ty < E) = E x 1 Ty<∞ exp - ∞ 0 L y T B y k(dy) > 0.
To prove (2.8), we fix t > 0. As above, a.s. on the event {t < τ c ∂ } = {σ t < T B 0 }, the function y → L y σt is continuous with compact support on (0, +∞). Since k is locally finite on (0, +∞), we deduce that κ t < ∞ a.s. Since in addition P x ( X t > 0) > 0 for all t > 0, (2.8) follows.

In the sequel, following the terminology of [START_REF] Feller | The parabolic differential equations and the associated semigroups of transformations[END_REF], we assume that +∞ is an entrance boundary, or equivalently (for diffusion processes on natural scale)

∞ 1 y(dm(y) + dk(y)) < ∞, (2.9)
and that 0 is either regular or an exit point, or equivalently

1 0 y(dm(y) + dk(y)) < ∞. (2.10)
Note that, in the case where k = 0, X = X is a diffusion without killing. In this case, the assumption ∞ 1 y dm(y) < ∞ corresponds to the fact that X comes down from infinity [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF], i.e. that there exist t > 0 and y > 0 such that inf x>y P x (T y < t) > 0, and

1 0 y dm(y) < ∞ is equivalent to assuming that τ ∂ = τ c
∂ < ∞ a.s. Further general properties of diffusion processes with killing are studied in Section 5, related to the notion of coming down from infinity.

Quasi-stationary distributions for diffusion processes with killing

We consider as above a diffusion process X on [0, +∞) with killing, on natural scale, with speed measure m(dx) and killing measure k(dx).

Exponential convergence to quasi-stationary distribution

We provide here sufficient criteria ensuring the existence of a unique quasistationary distribution for X, with exponential convergence of conditional distributions in total variation. Both criteria involve the diffusion process without killing X of Section 2.

Condition (C) Assume that

∞ 0 y (m(dy) + k(dy)) < ∞ and that there exist two constants t 1 , A > 0 such that

P x (t 1 < T 0 ) ≤ Ax, ∀x > 0 (3.1)
and

P x (τ d ∂ < τ c ∂ ) ≤ Ax, ∀x > 0. (3.2) Condition (C') Assume that ∞ 0 y (m(dy) + k(dy)
) < ∞ and that there exist three constants t 1 , A, ε > 0 such that

P x (t 1 < T 0 ) ≤ Ax, ∀x > 0
and k is absoluely continuous w.r.t. m on (0, ε) and satisfies

dk dm (x) ≤ A, ∀x ∈ (0, ε). (3.3)
We give practical criteria to check Conditions (C) and (C') is Subsection 3.2 and we give examples of applications and compare with existing results in Section 4.

Theorem 3.1. Assume that X is a one-dimensional diffusion on natural scale with speed measure m and killing measure k. If Assumption (C) or Assumption (C') is satisfied, then there exists a unique probability measure α on (0, +∞) and two constants C, γ > 0 such that, for all initial distribution µ on E,

P µ (X t ∈ • | t < τ ∂ ) -α(•) T V ≤ Ce -γt , ∀t ≥ 0. (3.4)
In this case, α is the unique quasi-stationary distribution for the process.

We recall that the quasi-stationary distribution α of Thm. 3.1 satisfies the following classical property: there exists λ 0 > 0 such that P α (t < τ ∂ ) = e -λ 0 t . The same assumptions also entail the following result on the asymptotic behaviour of absorption probabilities. Proposition 3.2. Assume that X is a one-dimensional diffusion on natural scale with speed measure m and killing measure k. If Assumption (C) or Assumption (C') is satisfied, then there exists a bounded function

η : (0, ∞) → (0, ∞) such that η(x) = lim t→∞ P x (t < τ ∂ ) P α (t < τ ∂ ) = lim t→+∞ e λ 0 t P x (t < τ ∂ ), (3.5)
where the convergence holds for the uniform norm. Moreover, α(η) = 1, η(x) ≤ Cx for all x ≥ 0 and some constant C > 0, and η belongs to the domain of the infinitesimal generator L of X on the set of bounded measurable functions on [0, +∞) equiped with uniform norm and

Lη = -λ 0 η. (3.6)
Still under the same assumptions, we also obtain the exponential ergodicity of the Q-process, which is defined as the process X conditioned to never be absorbed. Theorem 3.3. Assume that X is a one-dimensional diffusion on natural scale with speed measure m and killing measure k. If Assumption (C) or Assumption (C') is satisfied, then the three following properties hold true.

(i) Existence of the Q-process. There exists a family (Q x ) x>0 of probability measures on Ω defined by

lim t→+∞ P x (A | t < τ ∂ ) = Q x (A)
for all F s -measurable set A, for all s ≥ 0. The process (Ω, (F t ) t≥0 , (X t ) t≥0 , (Q x ) x>0 ) is a (0, ∞)-valued homogeneous strong Markov process.

(ii) Transition kernel. The transition kernel of the Markov process X under (Q x ) x>0 is given by p(x; t, dy) = e λ 0 t η(y) η(x) p(x; t, dy)1 y>0 , ∀t ≥ 0, x > 0, where p(x; t, dy) is the transition distribution of the diffusion with killing X.

(iii) Exponential ergodicity. The probability measure β on (0, +∞) defined by

β(dx) = η(x)α(dx).
is the unique invariant distribution of X under Q. Moreover, there exist positive constants C, γ such that, for any initial distributions µ 1 , µ 2 on (0, +∞),

Q µ 1 (X t ∈ •) -Q µ 2 (X t ∈ •) T V ≤ Ce -γt µ 1 -µ 2 T V ,
where

Q µ = (0,∞) Q x µ(dx).
The proofs of these results are given in Section 6.

On Conditions (C) and (C')

The first step to check Conditions (C) or (C') is to check (3.1). This condition is the one ensuring exponential convergence of conditional distributions of general diffusion processes without killing [4, Thm. 3.1], and it has been extensively studied in [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]. We recall here the most general criterion, which applies to all practical cases where ∞ is entrance and 0 is exit or natural, since it fails only when m has strong oscillations near 0 (see [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]Rk. 4]). Then, for all t > 0, there exists A t < ∞ such that

P x (t < T 0 ) ≤ A t x, ∀x > 0.
Using this result, the remaining conditions in (C) and (C') only deal with the behavior of the diffusion process X near 0. We provide two sufficient criteria for (3.2) in Condition (C) in Proposition 3.5 and Proposition 3.6. Note that checking property (3.3) in Condition (C') from given measures m and k is straightforward. Proposition 3.5. Assume that ∞ 0 y k(dy) < ∞ and that there exist ε > 0 and a constant C > 0 such that k is absolutely continuous w.r.t. the Lebesgue measure Λ on (0, ε) and such that

dk dΛ (x) ≤ C x , ∀x ∈ (0, ε). (3.8)
Then there exists a positive constant A > 0 such that

P x (τ d ∂ < τ c ∂ ) ≤ Ax.
Note that, contrary to Condition (C'), the previous criterion does not require the killing rate dk/dm of the diffusion process to be bounded on a neighbourhood of 0. It is actually independent of the measure m. For instance, if m(dx) = dx and k(dx) = 1

x dx on a neighbourhood of 0, then the assumptions of Proposition 3.5 are satisfied while dk/dm(x) = 1/x is unbounded near 0.

Proof of Proposition 3.5. We have from (2.3)

P x (τ d ∂ < τ c ∂ ) = P x κ T 0 > E = P x ∞ 0 L y σ T 0 dk(y) > E
Hence, by (2.2),

P x (τ d ∂ < τ c ∂ ) = P x ∞ 0 L y T B 0 dk(y) > E = P Y x T Y 0 > E ,
where Y is a diffusion process on (0, ∞) on natural scale with speed measure k, absorbed at 0 and without killing. We have

P Y x (T Y 0 > E) ≤ P Y x (E < T Y 0 < T Y ε ) + P Y x (T Y ε ≤ T Y 0 ) ≤ P Y x (E < T Y 0 < T Y ε ) + x/ε,
since Y is on natural scale and hence is a local martingale. As a consequence, it only remains to prove that

P Y x (E < T Y 0 < T Y ε ) ≤ Ax.
Because of our assumption (3.8),

P Y x (E < T Y 0 < T Y ε ) = P x T B 0 < T B ε ∩ ∞ 0 L y T B 0 dk(y) > E ≤ P x T B 0 < T B ε ∩ ∞ 0 L y T B 0 dy y > E/C ≤ P x ∞ 0 L y T B 0 dy y > E/C = 1 -v(x)
where v(x) = P Z x (T Z 0 ≤ E/C) with Z the diffusion on (0, +∞) on natural scale and with speed measure (1/y)dy, i.e. solution to the SDE

dZ t = Z t dW t ,
for (W t ) t≥0 a standard Brownian motion. The diffusion process (Z t ) t≥0 is a continuous state branching process (Feller diffusion), hence by the branching property,

v(x + y) = v(x)v(y), ∀x, y ≥ 0.
Since v(x) is non-increasing and v(0) = 1, we deduce from standard arguments that there exists a constant A > 0 such that v(x) = e -Ax . Proposition 3.5 follows.

Proposition 3.6. Assume that ∞ 1 yk(dy) < ∞ and k([0, 1]) < ∞. Then there exists a positive constant A > 0 such that

P x (τ d ∂ < τ c ∂ ) ≤ Ax.
We emphasize that measures k which are not absolutely continuous w.r.t. the Lebesgue measure are of course also covered by Prop. 3.6. Note also that Prop. 3.6 does not imply Prop. 3.5 since the killing measure dk(x) = (1/x)dΛ(x)1 x∈(0,1) is not finite.

Proof of Proposition 3.6. As in the beginning of the proof of Proposition 3.5, we obtain

P x (τ c ∂ < τ d ∂ ) = P x ∞ 0 L y T B 0 dk(y) < E = E x exp - ∞ 0 L y T B 0 dk(y) ≥ 1 -E x ∞ 0 L y T B 0 dk(y) = 1 -E Y x (T 0 | Y 0 = x),
where E Y denotes the expectation w.r.t. the law of a diffusion process Y on (0, +∞) on natural scale, without killing and whose speed measure is k(dy).

Finally, using classical Green formula for diffusion processes without killing applied to Y , we conclude that

P x (τ c ∂ < τ d ∂ ) ≥ 1- ∞ 0 (x∧y) k(dy) ≥ 1-xk([0, 1])-x ∞ 1 y k(dy) ≥ 1-Cx.
4 Examples and comparison with existing results

On general diffusions

Let us first recall that our results also cover the case of general (i.e. not necessarily on natural scale) killed diffusion processes (Y t , t ≥ 0) on [0, ∞)

with speed measure m Y and killing measure k Y , such that the diffusion without killing Y hits 0 in a.s. finite time. Under these assumptions, there exists a continuous and strictly increasing scale function s : [0, ∞) → [0, ∞) of the process Y such that s(0) = 0 and s(∞) = ∞, and our results apply to the process X t = s(Y t ) on natural scale, whose speed measure and killing measure are respectively given by In particular, we deduce from (4.1) the following criteria for (D) and (D').

m X (dx) = m Y * s(dx) and k X (dx) = k Y * s(dx), (4.1
P(t 1 < T 0 | Y 0 = y) ≤ As(y), ∀y > 0. Then (i) If k Y ([0, 1]) < ∞, then Y satisfies (D). (ii) If there exists ε, ρ, C > 0 such that dk Y dΛ (y) ≤ C s(y)
, ∀y ∈ (0, ε)

and s(y) -s(x) ≥ ρ(y -x), ∀0 < x < y < ε, (4.2)
then Y satisfies (D).

(iii) If the killing rate κ Y (y) := dk Y dm Y (y) exists and is uniformly bounded on (0, ε) for some ε > 0, then Y satisfies (D').

Proof. Point (i) follows immediately from Prop. 3.6 since k X ([0, s( 1

)]) = k Y ([0, 1]) < ∞.
For Point (ii), we deduce from (4.1) that

dk X d(Λ * s) ≤ C x , ∀x ∈ (0, s(ε)).
In view of of Prop. 3.5, it only remains to prove that d(Λ * s)/dΛ ≤ C in the neighborhood of 0. This is exactly (4.2). Finaly, under the conditions of Point (iii), by (4.1), κ X (y) := dk X dm X (x) exists and is uniformly bounded on (0, s(ε)). Hence X satisfies (C').

Let us also mention that, as will appear clearly in the proofs of Section 6, our methods can be easily extended to diffusion processes on a bounded interval, where one of the boundary point is an entrance boundary and the other is exit or regular, and also to cases where both boundary points are either exit or regular.

We provide now a few examples of diffusion processes that enter our setting (directly or considering the above extension).

Example 1. Let X be a sticky Brownian motion on the interval (-1, 1) which is sticky at 0, absorbed at -1 and 1, and with measurable killing rate κ : (-1, 1) → (0, +∞) such that (-1,1) κ(x)dx < ∞ and κ(0) < ∞. In this case X is the diffusion with killing on natural scale with speed measure and killing measure m(dx) = Λ(dx) + δ 0 (dx) and k(dx) = κ(x)Λ(dx) + κ(0)δ 0 (dx).

Extending our result to diffusion processes on (-1, 1) with regular boundaries, we immediately deduce from Proposition 3.6 that X admits a unique quasi-stationary distribution and satisfies (3.4).

Example 2. Let Y be a solution to the SDE

dY t = Y t dB t + Y t (r -cY t )dt,
where r ∈ R represents the individual growth rate without competition and c > 0 governs the competitive density dependence in the growth rate. We assume that the diffusion process Y is subject to the killing rate κ(y) = sin(1/y) ∨ √ y. Then Y admits a unique quasi-stationary distribution and satisfies (3.4) (see Proposition 4.2 for a detailed study of this case in a more general setting).

In this example, 0 is an exit boundary and κ is not continuous at 0 (see the comparison with the literature in Section 4.3).

Example 3. Let Y be a solution to the SDE

dY t = dB t -Y 2 t dt,
absorbed at 0. We may for example assume that the diffusion process Y is subject to the non-negative killing rate κ(y) = 1 √ y ∨ √ y, or to the killing measure

k(dy) = n∈N b n δ an (dy),
for any bounded sequence (a n ) n∈N and summable family (b n ) n∈N . In both cases, condition (D) is satisfied.

These two examples illustrate that our results also cover cases of unbounded (and even singular) killing rate κ.

Example 4. In addition, our method is general enough to apply to more complex processes. An example of diffusion without killing but with jumps is studied in [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]Sec. 3.5.4]. The extension to the case of diffusions with killing is straightforward using the same method: we consider a diffusion process on natural scale (X t , t ≥ 0) on [0, ∞) with speed measure m and killing measure k such that either Condition (C) or (C') is satisfied. Let us denote by L the infinitesimal generator of X.

Our first example is the Markov process ( X t , t ≥ 0) with infinitesimal generator

Lf (x) = Lf (x) + (f (x + 1) -f (x))1 x≥1 ,
for all f in the domain of L. In other words, we consider a càdlàg process following a diffusion process with speed measure m and killing measure k between jump times, which occur at the jump times of an independent Poisson process (N t , t ≥ 0) of rate 1, with jump size +1 if the process is above 1, and 0 otherwise. Then, a straightforward adaptation of Prop. 3.10 of [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF] implies that the conclusion of Thm. 3.1 holds for X.

Another simple example is given by the Markov process ( X t , t ≥ 0) with infinitesimal generator

Lf (x) = Lf (x) + f (x + 1) -f (x),
for all f in the domain of L. Here, +1 jumps occur at Poisson times, regardless of the position of the process. Then, in the case where (m, k) satisfies Condition (C) and if dm dΛ (x) ≤ C x in the neighborhood of 0, the conclusion of Thm. 3.1 holds for X.

To prove this, the only new difficulty compared with the previous case is to check that the probability that the process jumps before being killed or hitting zero is smaller than Cx when X 0 = x. By Condition (C), this is equivalent to prove the same property for the process without killing. Since the first jump time is exponential of parameter 1, this can be proved with the exact same argument as in Prop. 3.5 since

τ c ∂ = ∞ 0 L y T B 0 dm(y).
Of course, many easy extensions are possible, for example with jumps from x to g(x) at constant rate, where g is a non-decreasing function and might be 0 on some interval. In fact, even random jumps can be easily covered provided monotonicity properties of jump measures. This also covers situations with bounded non-decreasing jump rate. The case of general rates is more complicated but could also be attacked with our method.

Application to models of population dynamics

Quasi-stationary distributions for one-dimensional diffusion processes have attracted much interest as an application to population dynamics. In [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF][START_REF] Cattiaux | Competitive or weak cooperative stochastic Lotka-Volterra systems conditioned on non-extinction[END_REF], the authors consider the conditional behavior of diffusion processes conditioned not to hit 0. In [START_REF] Hening | A stochastic Lotka-Volterra Model with killing[END_REF], the authors consider one-dimensional diffusion processes that cannot hit 0 continuously (0 is a natural boundary) but are subject to a killing rate modeling the risk of catastrophic events. In this section, our aim is to bring together both points of view by considering logistic Feller diffusion processes and their extensions to general drifts, as in [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF], with an addition killing rate.

More precisely, let us consider a diffusion process Y on [0, +∞) solution to the SDE

dY t = Y t dB t + Y t h(Y t )dt,
where h is a measurable function from R + → R modeling the density dependence of the individual growth rate of the population. We assume that Y is subject to the nonnegative killing rate κ = dk Y dm Y ∈ L 1 loc ((0, +∞)). For such models, the conditions (D) and (D') directly give criteria for exponential convergence to the quasi-distribution. The next result gives explicit conditions on κ for a large class of functions h. Other conditions on κ could be obtained for other asymptotic behaviors of h. then there exists a unique probability measure α on (0, +∞) and two constants C, γ > 0 such that, for all initial distribution µ on (0, ∞),

P µ (X t ∈ • | t < τ ∂ ) -α(•) T V ≤ Ce -γt , ∀t ≥ 0.
Proof. In this case, the scale function and the speed measure associated with Y are given by s(y) = The killing measure associated to Y is k Y = κ(y)m(dy) = κ(y)e 2 y 0 h(u)du y dy. Our aim is to prove that condition (D) or (D') holds true using Proposition 4.1.

To do so, let us first check that

(0,+∞) s(y)(m Y (dy) + k Y (dy)) < ∞. (4.3)
We have

(0,+∞) s(y)(m Y (dy) + k Y (dy)) = (0,+∞) y 0 e -2 u 0 h(z)dz du (1 + κ(y))e 2 y 0 h(u)du y dy. Since h ∈ L 1 loc ([0, +∞)), we have y 0 e -2 u 0 h(z)dz du (1 + κ(y))e 2 y 0 h(u)du y ∼ y→0 1 + κ(y),
which is integrable in any bounded neighborhood of 0 by assumption. As a consequence,

(0,y 0 ) s(y)(m Y (dy) + k Y (dy)) < ∞.
On the other hand, e -2C/(1+β)(y 1+β -u 1+β ) du.

We have, for all y ≥ y 0 ,

y 0 0 e 2 y u h(z)dz du ≤ y 0 0 e 2 y 0 u h(z)dz du e -2C/(1+β)(y 1+β -y 1+β 0 ) = C e -2Cy 1+β /(1+β)
for some finite constant C . Moreover, setting

C β = sup x∈(0,1) (1 -x)/(1 - x 1+β ) < ∞, we have y y 0 e -2C/(1+β)(y 1+β -z 1+β )dz du ≤ y y 0 e -2C/(1+β)C -1 β (y-z)y β dz du ≤ (1 + β)C β 2Cy β .
As a consequence, there exists C > 0 such that

(y 0 ,+∞) s(y)(m Y (dy) + k Y (dy)) ≤ C (y 0 ,+∞) 1 + κ(y) y 1+β dy,
which is finite by assumption. Finally, we conclude that (4.3) holds.

The diffusion process X = s(Y ) is on natural scale with speed measure m X = m Y * s, which is dominated by (2/x) dΛ(x) in the neighborhood of 0, since s(y) ∼ 0+ y and s (y) → 0+ 1. We deduce from Thm. 3.4 that X satisfies (3.1) and hence that Y satisfies

P(t 1 < T 0 | Y 0 = y) ≤ As(y), ∀y > 0.
On the one hand, if

1 0 κ(y) y dy < ∞ then 1 0 k Y (dy) < ∞,
which allows us to conclude that (D) holds using Point (i) of Proposition 4.1. On the other hand, if lim sup x→0 κ(x) < ∞, then Point (iii) of the same proposition holds true and hence (D') holds.

On processes solutions to stochastic differential equations and comparison with the literature

In the case where the speed measure m is absolutely continuous w.r.t. the Lebesgue measure on (0, ∞), our diffusion processes on natural scale are solutions, before killing, to SDEs of the form

dX t = σ(X t )dB t , (4.4)
where σ is a measurable function from (0, ∞) to itself such that the speed measure m(dx) = 1 σ 2 (x) dx is locally finite on (0, ∞). Following the scale function trick of Section 4.1, our results actually cover all SDEs of the form

dY t = σ(Y t )dB t + b(Y t )dt such that b/σ 2 ∈ L 1 loc ((0, ∞)) (see Chapter 23 of [9]
). Existence of quasi-stationary distribution of diffusion processes with killing has been already studied in [START_REF] Steinsaltz | Quasistationary distributions for one-dimensional diffusions with killing[END_REF][START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]. These works are based on a careful study of the generator of the process and Sturm-Liouville methods. Their results cover the case of entrance or natural boundary at ∞, but only the case of regular boundary at 0. Our results only cover the case of entrance boundary at ∞, but 0 can be either natural or exit. In the case where ∞ is an entrance boundary, the stronger previously known result is [START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]Thm. 4.13] for diffusion processes (Y t , t ≥ 0) on [0, ∞) which are solution to SDEs of the form

dY t = dB t -h(Y t )dt, (4.5)
with killing measure absolutely continuous w.r.t. the Lebesgue measure with continuous density on [0, ∞). Let us recall that the scale function of the solution (Y t , t ≥ 0) to (4.5) is given by

s(x) = x 0 exp 2 y 1 h(z) dz dy,
and that X t = s(Y t ) is solution to the SDE

dX t = exp 2 s -1 (Xt) 1 h(z) dz dB t . (4.6)
In particular, this diffusion with killing is not as general as (4.4) since it only covers the case where the diffusion coefficient is positive and continuous. In fact, in all the references previously cited, the drift h is assumed at least continuous so that the diffusion coefficient in (4.6) is at least C 1 . In addition, our criteria allow unbounded killing rates at 0 as shown in the next proposition, contrary to existing results. This result can be proved similarly as Proposition 4.2. In the case where κ is continuous bounded and h is continuous on [0, +∞), existence of a quasi-stationary distribution was already proved in [START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF]. However, even in this case, uniqueness of the QSD, exponential convergence of conditional distributions and exponential ergodicity of the Q-process are new. Conversely, other results of [START_REF] Kolb | Quasilimiting behavior for one-dimensional diffusions with killing[END_REF] cannot be easily studied with our methods, such as dichotomy results on the behavior of conditioned diffusion processes with killing and the case of natural boundary at ∞.

On continuous and discontinuous absorption times

when 0 is regular or exit and ∞ is entrance

We give in this section a result on the probability P x (τ c ∂ < t ∧ τ d ∂ ) that the process hits 0 continuously before being killed, which is useful for the proof of our results of Subsection 3.1. In the introduction, we called this property coming down from infinity before killing in reference to the classical property of coming down from infinity for diffusion processes without killing (see [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF][START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]). Since this result also has an interest by itself, we give it in a separate section.

Theorem 5.1.

∃t > 0 such that lim x→+∞ P x (τ c ∂ < t ∧ τ d ∂ ) > 0 ⇐⇒ ∞ 0 y (dk(y) + dm(y)) < ∞. (5.1)
This result is obtained as a consequence of the next proposition.

Proposition 5.2. For all k and m and all t > 0,

P x (τ c ∂ < t ∧ τ d ∂ ) = P x (τ c ∂ < t) = E 1 ∞ 0 Z x s m(ds)<t exp - ∞ 0 Z x s k(ds) , (5.
2) where

Z x t = 2 t 0 Z x s dW s + 2(t ∧ x), ∀t ≥ 0.
In particular, P x (τ c ∂ < t ∧ τ d ∂ ) is non-increasing w.r.t. x > 0 and

lim x→+∞ P x (τ c ∂ < t ∧ τ d ∂ ) = E 1 ∞ 0 Zs m(ds)<t exp - ∞ 0 Z s k(ds) , (5.3)
where (Z t , t ≥ 0) is a squared Bessel process of dimension 2 started from 0.

Proof. We have

P x (τ c ∂ < t ∧ τ d ∂ ) = P x ( T 0 < t and κ T 0 < E) = P x ∞ 0 L y σ T 0 m(dy) < t and ∞ 0 L y σ T 0 k(dy) < E = P x ∞ 0 L y T B 0 m(dy) < t and ∞ 0 L y T B 0 k(dy) < E = E x 1 ∞ 0 L y T B 0 m(dy)<t exp - ∞ 0 L y T B 0 k(dy) .
(A2) there exists c 2 > 0 such that for all x > 0 and t ≥ 0,

P ν (t < τ ∂ ) ≥ c 2 P x (t < τ ∂ ).
Hence, we need first to prove that (C) implies (A) (in Subsection 6.1), second, to prove that (C') implies (A) (in Subsection 6.2), and finally, to prove that η(x) ≤ Cx for all x ≥ 0 (in Subsection 6.3), which is the only part of Proposition 3.2 left to prove.

Proof that Condition (C) implies (A)

Note that [3, Thm. 2.1] also assumes that

P x (t < τ ∂ ) > 0, ∀x > 0, ∀t > 0, (6.1)
which is entailed by Prop. 2.2. Our proof of (A) follows four steps, similarly as the proof of [START_REF] Champagnat | Exponential convergence to quasistationary distribution for one-dimensional diffusions[END_REF]Thm. 3.1]. In the first step, we prove that when X 0 is close to 0, then, conditionally on non-absorption, the process exits some neighborhood of 0 with positive probability in bounded time. In the second step, we construct the measure ν involved in (A). We next prove (A1) in the third step, and (A2) in the last step.

Step 1: the conditioned process escapes a neighborhood of 0 in finite time. The goal of this step is to prove that there exists ε, c > 0 such that

P x (X t 1 ≥ ε | t 1 < τ ∂ ) ≥ c, ∀x > 0, (6.2)
where t 1 is taken from Assumption (C).

To prove this, recall that X is a diffusion process on natural scale with speed measure m but with null killing measure. We first observe that, since X is a local martingale, for all x ∈ (0, 1),

x = E x ( X t 1 ∧ T 1 ∧τκ ),
where we recall that τ

κ = τ d ∂ on the event τ d ∂ < ∞ and T 0 = τ c ∂ on the event τ c ∂ < ∞.
But the absorption at 0 ensures that

X t 1 ∧ T 1 ∧τκ = X t 1 ∧ T 1 ∧τκ 1 t 1 ∧τκ< T 0 + 1 T 1 < T 0 <t 1 ∧τκ .
We thus have for all x ∈ (0, 1)

x = P x (t 1 ∧ τ κ < T 0 )E x ( X t 1 ∧ T 1 ∧τκ | t 1 ∧ τ κ < T 0 ) + P x ( T 1 < T 0 < t 1 ∧ τ κ ).
Now, the strong Markov property entails

P x ( T 1 < T 0 < t 1 ∧ τ κ ) ≤ P x ( T 1 < T 0 )P 1 ( T 0 < t 1 ∧ τ κ ).
Since P 1 ( T 0 < t 1 ∧ τ κ ) < 1 (see Prop. 5.2) and P x ( T 1 ≤ T 0 ) = x (recall that X is a local martingale), we deduce from the two previous equations that there exists a constant A > 0 such that

A x ≤ P x (t 1 ∧ τ κ < T 0 )E x ( X t 1 ∧ T 1 ∧τκ | t 1 ∧ τ κ < T 0 ).
By Assumption (C), we have

P x (t 1 < T 0 ) ≤ Ax and P x (τ d ∂ < τ c ∂ ) ≤ Ax. But, by definition of τ c ∂ and τ d ∂ , we have {τ d ∂ < τ c ∂ } = {τ κ < T 0 }, so that P x (t 1 ∧ τ κ < T 0 ) ≤ 2Ax.
As a consequence,

E x 1 -X t 1 ∧ T 1 ∧τκ t 1 ∧ τ κ < T 0 ≤ 1 - A 2A .

Since we can assume without loss of generality that

A 2A < 1, Markov's in- equality implies that, setting b = 1 -1 -A /2A, P x X t 1 ∧ T 1 ∧τκ ≤ b t 1 ∧ τ κ < T 0 ≤ 1 -A /2A 1 -b = 1 -b, hence P x X t 1 ∧ T 1 ∧τκ > b t 1 ∧ τ κ < T 0 ≥ b.
We deduce that, for all x ∈ (0, b),

P x ( T b < t 1 ∧ τ κ | t 1 ∧ τ κ < T 0 ) = P x ( T b < t 1 ∧ T 1 ∧ τ κ | t 1 ∧ τ κ < T 0 ) ≥ b. (6.3) Now, since P b (t 1 < T 0 ∧ τ κ ) > 0 (see Prop. 2.2), there exists ε ∈ (0, b) such that P b (t 1 < T ε ∧ τ κ ) > 0. (6.4)
Hence, we deduce from the strong Markov property that

P x (X t 1 ≥ ε) = P x ( X t 1 ≥ ε and t 1 < τ κ ) ≥ P x ( T b < t 1 ∧ τ κ )P b (t 1 < T ε ∧ τ κ ) ≥ b P x (t 1 ∧ τ κ < T 0 ) P b (t 1 < T ε ∧ τ κ ),
where the last inequality follows from (6.3). But P

x (t 1 ∧ τ κ < T 0 ) ≥ P x (t 1 < T 0 ∧ τ κ ) = P x (t 1 < τ ∂ ), so that P x (X t 1 ≥ ε | t 1 < τ ∂ ) ≥ b P b (t 1 < T ε ∧ τ κ ) > 0.
This entails (6.2) for x < b. For x ≥ b, the continuity (before killing) and the strong Markov property for X imply

P x (X t 1 > ε | t 1 < τ ∂ ) ≥ P x (X t 1 > ε) ≥ P x (T b < ∞)P b (t 1 < T ε ∧ τ κ ) ≥ P x (τ c ∂ < ∞)P b (t 1 < T ε ∧ τ κ )
Using Proposition 5.2 and (6.4), the proof of (6.2) is completed.

Step 2: Construction of coupling measures for the unconditioned process.

We prove that there exist two constants t 2 , c > 0 such that, for all x ≥ ε,

P x (X t 2 ∈ •) ≥ c ν, (6.5)
where ν

= P ε (X t 2 ∈ • | t 2 < τ ∂ ).
Fix x ≥ ε and construct two independent diffusions X ε and X x with speed measure m(dx), killing measure k(dx), and initial value ε and x respectively. Let θ = inf{t ≥ 0 : X ε t = X x t }. By the strong Markov property, the process

Y x t = X x t if t ≤ θ, X ε t if t > θ
has the same law as X x . By the continuity of the paths of the diffusions before killing, θ ∧ τ d,x ∂ ≤ τ x ∂ := inf{t ≥ 0 : X x t = 0}. Hence, for all t > 0,

P(θ < t) ≥ P(τ c,x ∂ < t and τ c,x ∂ < τ d,x ∂ ).
By Theorem 5.1, there exists t 2 > 0 and c > 0 such that

inf y>0 P y (τ c ∂ < t 2 and τ c ∂ < τ d ∂ ) ≥ c > 0. Hence P x (X t 2 ∈ •) = P(Y x t 2 ∈ •) ≥ P(X ε t 2 ∈ •, τ c,x ∂ < t 2 and τ c,x ∂ < τ d,x ∂ ) ≥ c P ε (X t 2 ∈ •),
where the last inequality follows from the independence of X x and X ε . Therefore, (6.5) is proved with c = c P ε (t 2 < τ ∂ ).

Step 3: Proof of (A1). Using successively the Markov property, Step 2 and Step 1, we have for all x > 0

P x (X t 1 +t 2 ∈ • | t 1 + t 2 < τ ∂ ) ≥ P x (X t 1 +t 2 ∈ • | t 1 < τ ∂ ) ≥ ∞ ε P y (X t 2 ∈ •)P x (X t 1 ∈ dy | t 1 < τ ∂ ) ≥ c ν(•)P x (X t 1 ≥ ε | t 1 < τ ∂ ) ≥ c c ν(•).
This entails (A1) with t 0 = t 1 + t 2 and c 1 = cc .

Step 4: Proof of (A2). Let a > 0 be such that ν([a, +∞)) > 0. Then we have, for all x ≥ a, Using this inequality and a standard renewal argument, we deduce that there exists ρ > 0 such that, for all k ∈ N, P a (X kt 0 ≥ a) ≥ e -ρkt 0 . (6.6)

Now, we know from Step 4 of the proof of [4, Theorem 3.1] (see also [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF]) that, since the diffusion X on [0, ∞) without killing has ∞ as an entrance boundary, for all r > 0, there exists y r > 0 (which can be assumed larger than a without loss of generality) such that sup x≥yr E x (e r Ty r ) < +∞.

This implies sup

x≥yρ E x (e ρ(Ty ρ ∧τ ∂ ) ) < +∞, (6.7)

for the constant ρ of ( Now, for all x ≥ y ρ , with a constant C > 0 that may change from line to line, using successively (6.7), (6.8), (6.9) and (6.6), we obtain P x (t < τ ∂ ) ≤ P x (t < T yρ ∧ τ ∂ ) + ≤ Ce -ρ t/t 0 t 0 + CP a (t < τ ∂ ) Finally, since ν([a, +∞)) > 0, (6.8), (6.10) and (6.11) entail (A2).

check that there exists a constant C such that η(x) ≤ Cx for all x ≥ 0, where η(x) is the uniform limit of e λ 0 t P x (t < τ ∂ ). Since for all t ≥ t 1 , where the constant t 1 is the one given in Assumption (C) or (C'), e λ 0 t P x (t < τ ∂ ) = e λ 0 t E x P Xt 1 (t -t 1 < τ ∂ ) ≤ Ce λ 0 t 1 P x (t 1 < τ ∂ ). (6.14) For the last inequality, we used the fact that e λ 0 s P y (s < τ ∂ ) is uniformly bounded in y > 0 and s ≥ 0 because of the uniform convergence in (3.5). By Assumption (C) or (C'), P x (t 1 < τ ∂ ) ≤ P x (t 1 < T 0 ) ≤ Ax, and we obtain the inequality η(x) ≤ Cx by letting t → +∞ in (6.14).

Theorem 3 . 4 ([ 4 ]

 344 Thm 3.7). Assume that ∞ 0 y m(dy) < ∞ and y) z 2 m(dz) dx < ∞.(3.7)

  ) where m Y * s and k Y * s denote the pushforward measures of m Y and k Y through the function s (the formula for m X can be found in [17, Thm. VII.3.6] and for the measure k X , one may use for example (2.3) and the fact that L y, Y t = L s(y), X t which follows from the occuptation time formula of [6, p. 160]). Hence X satisfies (C) or (C') if and only if Y satisfies the following Conditions (D) or (D') respectively. In particular, the conclusions of Theorem 3.1 apply to Y if (D) or (D') is fulfilled. Condition (D) Assume that ∞ 0 s(y) (m Y (dy) + k Y (dy)) < ∞. Assume also that there exist two constants t 1 , A > 0 such that P(t 1 < T 0 | Y 0 = y) ≤ As(y), ∀y > 0 and P(τ d ∂ < τ c ∂ | Y 0 = y) ≤ As(y), ∀y > 0, where the absorption times T 0 , τ d ∂ and τ c ∂ are constructed here from the processes Y and Y . Condition (D') Assume that ∞ 0 s(y) (m Y (dy) + k Y (dy)) < ∞. Assume also that there exist three constants t 1 , A, ε > 0 such that P(t 1 < T 0 | Y 0 = y) ≤ As(y), ∀y > 0 and k Y is absoluely continuous w.r.t. m Y on (0, ε) and satisfies dk Y dm Y (y) ≤ A, ∀y ∈ (0, ε).

Proposition 4 . 1 .

 41 Assume that ∞ 0 s(y) (m Y (dy) + k Y (dy)) < ∞ and that there exist two constants t 1 , A > 0 such that

Proposition 4 . 2 .

 42 Assume that h ∈ L 1 loc ([0, +∞)) and that there exists C, β > 0 such that, for some y 0 > 0, h(y) ≤ -Cy β , ∀y ≥ y 0 .

y 0 e

 0 -2 u 0 h(z)dz du and m Y (dy) = e 2 y 0 h(u)du y dy.

Proposition 4 . 3 .

 43 Assume that h ∈ L 1 loc ([0, +∞)) and that there exists C > 0 and β > 1 such that, for some y 0 > 0, h(y) ≤ -Cy β , ∀y ≥ y 0 . y β dy < ∞, then (D) is satisfied.

P

  x (X t 0 ≥ a) ≥ c 1 ν([a, +∞))P x (τ ∂ > t 0 ) ≥ c 1 ν([a, +∞))P x (T a < ∞)P a (τ ∂ > t 0 ) ≥ c 1 ν([a, +∞))P x (τ c ∂ < τ d ∂ )P a (τ ∂ > t 0 ),where we have used the strong Markov property for the second inequality. By Theorem 5.1, inf x∈(0,+∞) P x (τ c ∂ < τ d ∂ ) > 0 and we know that P a (τ ∂ > t 0 ) > 0 by Proposition 2.2. As a consequence, inf x≥a P x (X t 0 ≥ a) > 0.

t 0 P

 0 yρ (t -s < τ ∂ ) P x (T yρ ∈ ds) ≤ Ce -ρt + C t 0 P a (t -s < τ ∂ ) P x (T yρ ∈ ds)

  6.6). Using the strong Markov property, for all x ∈ [a, +∞) and y ∈ [a, y ρ ],P x (t < τ ∂ ) ≥ P x (T a < ∞)P a (T y < ∞)P y (t < τ ∂ ) ≥ P x (τ c ∂ < τ d ∂ )P a (T yρ < ∞)P y (t < τ ∂ ).Using Theorem 5.1 and Proposition 2.2, we infer that there exists C > 0 such that, for all t ≥ 0, sup x∈[a,yρ]P x (t < τ ∂ ) ≤ C inf x∈[a,+∞) P x (t < τ ∂ ). (6.8)Finally, we also deduce from the Markov property that, for all s < t, P a (X s/t 0 t 0 ≥ a) inf

x∈[a,+∞)

P x (t -s < τ ∂ ) ≤ P a (t < τ ∂ ).

(6.9)

  (X s/t 0 t 0 ≥ a) P Now, for all x ∈ [a, +∞) and all y ∈ (0, a),P x (t < τ ∂ ) ≥ P x (T y < ∞)P y (t < τ ∂ ) ≥ P x (τ c ∂ < τ d ∂ )P y (t < τ ∂ ).

	Using Proposition 5.2, we deduce that there exists C > 0 such that, ∀t ≥ 0,
	sup x∈(0,a)	P x (t < τ ∂ ) ≤ C inf x∈[a,+∞)	P x (t < τ ∂ ).	(6.11)

t 0 1 P a x (T yρ ∈ ds) ≤ CP a (t < τ ∂ ) + CP a (t < τ ∂ ) t 0 e ρs P x (T yρ ∈ ds).

We deduce that, for all t ≥ 0, sup x∈[yρ,+∞)

P x (t < τ ∂ ) ≤ C inf x∈[a,+∞) P x (t < τ ∂ ).

(6.10)

Hence (5.2) follows from Ray-Knight's theorem (cf. e.g. [START_REF] Rogers | Diffusions, Markov processes, and martingales[END_REF]).

Since the process Y x t = Z x t is solution to Proof of Theorem 5.1. In view of (5.3), the equivalence (5.1) follows from [START_REF] Pitman | A decomposition of Bessel bridges[END_REF]Prop. (2.2), Lemma (2.3)], which states that, for any square Bessel process (R t , t ≥ 0) of positive dimension started from 0 and any positive Radon measure µ on (0, +∞),

Hence both r.v. Assumption (A) There exists a probability measure ν on (0, +∞) such that (A1) there exists t 0 , c 1 > 0 such that for all x > 0,

Proof that Condition (C') implies (A)

The proof follows exactly the same steps as in Subsection 6.1. The only step which needs to be modified is Step 1, so we only detail this step.

Step 1: the conditioned process escapes a neighborhood of 0 in finite time. The goal of this step is to prove that there exists ε , c > 0 such that

(6.12)

We know from the study of diffusions without killing of [4, Step 1 of Section 5.1] that there exists c, b > 0 such that, for all x ≤ b,

By Assumption (C') and using (2.1), a.s. for all t ≤ T ε ,

Note that we can assume without loss of generality in (6.13) that b ≤ ε.

Therefore, for all x ∈ (0, b), a.s. under P x ,

Hence, using (6.13) and the independence between E and X,

But, for all ε ∈ (0, b), x ∈ (0, b),

By Proposition 2.2, the last factor of the r.h.s. is positive for ε > 0 small enough, which concludes the proof of (6.12), since the case x ≥ b can be handled exactly as in Subsection 6.1.

6.3 Proof that η(x) ≤ Cx for all x ≥ 0

We recall that, since we already proved that (C) and (C') imply (A), [3, Prop. 2.3] entails the main part of Proposition 3.2, and we only have to


Fabien Tarissan

Leo Liberti

Camilo La Rota

Network reconstruction: a mathematical programming approach

We propose a method employing mathematical programming and global optimization techniques for solving inverse problems arising in biological regulatory network (BRN) reconstruction. This problem consists in estimating unknown parameters of a model that describe the structure and dynamics of a biological system from a set of experimental observations and can be naturally cast as an optimization problem: choose the parameter values minimizing a given distance between the observed and estimated values of some observable variables. This minimization is subject to constraints derived from the models.

Two significant examples are presented on how to handle different kinds of dynamics: pattern formation occurring from diffusable and nondiffusable gene products in the drosophila melanogaster morphogenesis and reconstruction of the gene regulatory network of arabidospsis thaliana based on the identification of stable sub-networks during morphogenesis.

Introduction

A typical problem in biological systems consists in the simulation of a biological process, where from certain initial or boundary conditions one proceeds to solve a system of equations describing the biological process and obtains trajectories describing the system evolution in time and/or space. This is what is known as a forward problem. An example of a forward problem is the computation of the time evolution of a Gene Regulatory Network (GRN) configuration (values associated to the genes and representing gene products concentrations) given an initial condition and some local interaction rules. The GRN may be described in a detailed-quantitative manner by a system of Nonlinear Differential-Algebraic Equations (NDAE), where the values associated to the genes are real and directly represent the gene products concentrations, or in a coarse-qualitative manner by a system of nonlinear Binary or Multivalued Difference Equations (NBdE,NMdE) or by a system of Boolean or Multivalued Logical Equations (LE,MLE), where the set of values represent discrete levels of gene products concentrations such as abscence (0) or saturation [START_REF] Balanzá | Patterning the female side of arabidopsis: the importance of hormones[END_REF] in the case of boolean representations.

1 That, however, presupposes that all "known values" (coefficients of the DAE system, initial concentration of gene products) are known apriori, which is rarely the case. Thus, pre-simulation methods addressing the problem of estimating the unknown parameters from a set of experimental observations become of paramount importance. Such problems are often termed inverse problems as in some sense one needs to solve a problem having the same form, but where the roles of known and unknown values are interchanged and because only the model equations for the forward problem are known.

The usual approach developped in the context of molecular biology in order to adress the problem is often based on simulation techniques. It requires first to design specific heuristics that will infer values to the missing parameters and then, in a second step, to simulate the evolution of the system in order to compare the values generated by the simulation with the observed ones. The gap between generated and observed values then defines the quality of the estimation. As one might guess, the efficiency of this approach is directly related to the heuristic used which induces generally to spend time and effort on tuning the algorithms for each study. This remark has a great impact as soon as we intend to cover a wide range of organisms.

The approach we propose in this paper relies on the remark that the problem of estimating unknown parameters from a set of experimental observations can be naturally cast as an optimization problem: choose the parameter values minimizing the sum of a given norm of the differences between the observed and estimated values ofsome observables. Although designing specific algorithms for a given organism could provide faster results, the approach we propose has the advantage of generality: potentially any parameter estimation problem can be modeled by mathematical programming and solved with corresponding generalpurpose algorithms. More precisely, we use the modeling language AMPL [START_REF] Fourer | The AMPL Book[END_REF] in order to describe the biological system and then apply different reformulations (both exact and approximative) in order to ease the computation process carried out by dedicated Non Linear Programs (NLP) and Mixed Integer NLP (MINLP) numerical solvers [START_REF] Liberti | Writing global optimization software[END_REF]. The relevance of a general-purpose approach that sacrifices some computational efficiency is to provide general modeling environments and software packages with solution methods that can easily adapt to all (or at least most) models input by the user. This need is in our case explicitly established by the EU Morphex project and its associated modeling language.

The rest of the paper is organized as follow. In Section 2 we present the basic concepts underlying mathematical programming. In Section 3 we show how to use such a technique to reconstruct the GRN of Drosophila meloganaster in the early development stage relying on NDAE. In Section 4 we test the same approach but on the Arabidopsis thaliana modeled using BNdE. Finally, in Section 5 we discuss the main advantages and drawbacks of our approach.

Mathematical programming and optimization

A mathematical programming problem is formulated as follows:

min x f (x) subject to g(x) ≤ 0, (1)
where x ∈ R n are the decision variables and f : R n → R is the objective function to be minimized subject to a set of constraints g : R n → R m which may also include variable ranges or integrality constraints on the variables. Global Optimization is concerned with the solution of problems (1) where f, g are nonconvex nonlinear forms. A problem where f, g are nonlinear is known as a Nonlinear Programming problem (NLP); if some integrality constraints are present on the variable bounds the problem is known as a Mixed-Integer NLP (MINLP).

A mathematical programming problem which has some integer variables but whose objective function and constraints are linear forms is called a Mixed-Integer Linear Programming (MILP) problem. There exists general-purpose solution algorithms, both exact and heuristic, for all problem forms in NLP, MINLP, MILP. Currently, MILP solution methods are the most advanced, and the de facto standard solver is the ILOG CPLEX [START_REF] Ilog | ILOG CPLEX 8.0 User's Manual[END_REF] solver. MINLPs and nonconvex NLPs can be solved by many different global optimization methods such as BARON [START_REF] Sahinidis | BARON 7.2.5: Global Optimization of Mixed-Integer Nonlinear Programs[END_REF]. Most frequently, NLPs and MINLPs undergo a reformulation stage before being solved [START_REF] Liberti | Writing global optimization software[END_REF][START_REF] Liberti | Reformulation techniques in mathematical programming[END_REF].

Reproducing a continuous regulation of genes

The motivation of the present section is to test the relevancy of this framework on a first biological case study. We focused here on the Gap-Gene formation in the early development of Drosophila melanogaster. At this stage (between cycle 13 and 14A), the embryo contains static nuclei which can diffuse their gene products within their neighborhood. The model [START_REF] Jaeger | Dynamical analysis of regulatory interactions in the gap gene system of drosophila melanogaster[END_REF][START_REF] Nanfack | Efficient parameter estimation for models of pattern formation in early embryogenesis of Drosophila melanogaster[END_REF] uses a discrete representation of the space (linear topology) but a continuous regulation of the genes in time. The differential equation that describe this continuous regulation is as follows:

dg ia (t) dt = R a Φ(u ia (t)) -λ a g ia (t) + D a (g i+1,a (t) -2g ia (t) + g i-1,a (t)), (2)
where g ia (t) is the concentration of gene a (belonging to the set of genes N γ ) in nucleus i (belonging to the set of nuclei N ι ) at time t, R a is the production rate for gene a, Φ(y) is the sigmoid regulation function

Φ(y) = 1 2 y y 2 + 1 + 1 (3)
defined for all y, u ia (t) is a function defined as follows:

u ia (t) = b∈N γ W ba g ib (t) + m a g bcd i + h a (4)
where W ba is the weight on the arc (b, a) in the digraph representation G = (N γ , A) of the gene regulatory network, m a is the regulatory influence of the maternal gene bcd, h a is the activation threshold for Φ; λ a is the decay rate and D a is the diffusion coefficient for gene a.

The solution method proposed in [START_REF] Nanfack | Efficient parameter estimation for models of pattern formation in early embryogenesis of Drosophila melanogaster[END_REF] is a Global Optimization (GO) evolutionary algorithm applied to the following unconstrained optimization problem:

min i∈N ι t (g ia (t) -g data ia (t)) 2 + Π R + Π λ + Π D + Π u , (5)
where g data ia (t) is a set of experimentally observed values for the concentration level of gene a in nucleus i at time t; Π R is a penalty function that helps restricting R a to lie in a pre-determined range [R L a , R U a ] (Π λ and Π D are similar to Π R ); and finally, Π u is defined as

Π u = e Θ -1 (6)
(where

Θ = Λ( (b,a)∈A (W ba v max b ) 2 + (m a v max bcd ) 2 + h 2 a )
) if Θ > 0 and 0 otherwise, and v max b and v max bcd are the maximum values for gene b ∈ N γ and bcd respectively. We remark straight away that since Θ is a sum of squares, it is necessarily nonnegative, which implies that it suffices to describe Π u simply as e Θ -1 without any need for conditional definitions. Consider now the functions α 1 (y) = y + 1 and α 2 (y) = log(y) for all y > 0. Since both are monotonically increasing, the minimization of Π u determines the same set of optimal solutions as the minimization of α 2 (α 1 (Π u )). Hence ( 5) is equivalent to:

min i∈N ι t (g ia (t) -g data ia (t)) 2 + Π R + Π λ + Π D + Θ. (7)
Thus, the model used in this example contains two different kinds of elements: the description of the equations (rules 2, 3, 4) that control the evolution of the parameters and the description of the criterion that allows to establish if the estimated values are close the real ones (rule 7). This distinction will be naturally reflected in the model described in AMPL. The rule 5 is for instance clearly the objective function we are looking for whereas all the equations that describes the dynamic of the model will be stated as constraints. In order to ease the reading, we don't present the complete model in AMPL in this section but in the appendix 6.1 . The only main difference between the two versions that worth being noticed is related to the way we model the terms Π R , Π λ and Π D . Since those functions corresponds to specific boundaries for the parameters R, λ and D, it is natural to use constraints as follow:

∀a ∈ N γ    R L ≤ R a ≤ R U λ L ≤ λ a ≤ λ U D L ≤ D a ≤ D U
It induces in particular that the objective function will have a form slightly different from the equation [START_REF] Nanfack | Efficient parameter estimation for models of pattern formation in early embryogenesis of Drosophila melanogaster[END_REF]:

min a∈N γ i∈N ι t∈T data (g a i (t) -g data a i (t)) 2 + a∈N γ b∈N γ (W a b v b max ) 2 + a∈N γ ((m a v bcd max ) 2 + h 2 a ). ( 8
)
4 Using fixed points for GRN reconstruction

In the previous section, we showed an example where the objective function derived from the equation ( 7) used to establish the quality of the estimated values can easly be translated into AMPL [START_REF] Sahinidis | BARON 7.2.5: Global Optimization of Mixed-Integer Nonlinear Programs[END_REF]. This is of course directly related to the way the model is presented. In order to validate our approach, it is necessary to experience the same process on another kind of criterion. This is what this section is devoted to. More precisely, we solve here an inverse problem consisting in the determination of the parameters of a GRN in such a way that the reconstructed network show observed states at given development stages. Those transitory states of the GRN are expressed as stable states of sub-networks such that the activities of the nodes involved in a sub-network do not change anymore as soon as the stable state is reached and the subnetwork is isolated from the rest of the network. The case of study chosen here is the Arabidopsis thaliana floral meristem GRN and the expression patterns observed at the very first development stages. This network is still poorly understood, even if data concerning genes and morphogenes involved, single regulatory interactions and small circuits exist in the literature [START_REF] Balanzá | Patterning the female side of arabidopsis: the importance of hormones[END_REF].

Given a directed graph G = (V, A), a discrete set of time instants T , a set of development stages S (both of which we suppose to be an initial contiguous proper subset of N) and the following functions:

• a function α : A → {+1, -1} called the arc sign function;

• a function ω : A → N called the arc weight function;

• a function x : V × T → {0, 1} called the gene activation function;

• a function ι : V → {0, 1} called the initial configuration;

• a function θ : V → Z called the threshold function, A gene regulatory network (GRN) is a 7-tuple (G, T, α, ω, x, ι, θ) such that:

∀v ∈ V x(v, 1) = ι(v) (9) ∀v ∈ V, t ∈ T {1} x(v, t) = 1 if u∈δ -(v) α(u, v)ω(u, v)x(u, t -1) ≥ θ(v) 0 otherwise, ( 10
)
where δ -(v) = {u ∈ V | (u, v) ∈ A} for all v ∈ V . Eqns. ( 9)-(10) together are called the evolution rules of the GRN. For any particular t ∈ T , x(•, t) : V → {0, 1} is called a configuration. Since the evolution rules relate a configuration at time t with a configuration at time t -1, if x(•, t) = x(•, t -1) then x(•, t ′ ) = x(•, t) for all t ′ > t: such configurations are called fixed points of the GRN.

In this section we address the following problem:

Stable Subnetworks Reconstruction of a GRNs (SSRGRN). Given a digraph G = (V, A), a time instant set T , a set S of development stages, a set R s of observed cellular types at each stage, an arc sign function α, a set of initial configurations for each stage and for each observed cellular type I rs : V → {0, 1}, a set U s ⊆ V determining the nodes of the (induced) subnetwork G s at each stage to be reconstructed (U s ⊆ U s+1 ⊆ V ) and the observed configurations O rs : U s → {0, 1}, find an arc weight function ω and a threshold function θ with the property that for all s ∈ S and all ι ∈ I rs there exists a gene activation function x such that (G s , T, α, ω, x, ι, θ) is a GRN subnetwork whose fixed points are at a minimum distance to the observed data O rs .

In other words, we attempt to estimate the arc weights and threshold functions of a GRN from the knowledge of the digraph topology G, the induced subnetwork topology U and the arc sign function α in such a way that (a) the GRN evolution rules are consistent with respect to a certain set of initial configurations and (b) the fixed points in the subnetwork induced by the estimated values are as close as possible to the given ones.

Our primary concern in solving the SSRGRN is thus modellistic rather than algorithmic. One of the foremost difficulties is that of employing a static modelling paradigm -such as mathematical programming -in order to describe a problem whose very definition depends on time. Another important difficulty resides in expressing the necessary and sufficient conditions for a configuration to be a fixed point as constraints suitable for use in a formulation like [START_REF] Balanzá | Patterning the female side of arabidopsis: the importance of hormones[END_REF]. In appendix 6.2, we present a simplified version in which only one targeted configuration S is used.

Conclusion

As seen in this paper, the approach we propose present the advantage of presenting a unified framework able to handle two different kinds of organisms. The chosen organisms are such that both the dynamics induced by the networks (continuous in the case of drosophila and discrete in arabidopsis) and the criterion used to drive the way the values are infered tends to prove that this approach covers any parameter estimation problem as soon as the evolution rules can be described as constraints. This provides therefore an attractive environment for solving inverse problems in general.

However, the arabidopsis case showed that the definition of the objective function must be delt carefully in order to render properly the dynamics that one wants to capture. It is then compelling to try to automatize the way we define the objective function given a set of crieria (mimicking continuous regulation, identifying sub-networks which are stable, etc. . . ). We are already engaged in this direction.

Besides, the description in AMPL of the model might also rise problems for the solvers (the definition of Φ in (3 for instance) and it is then tempting to automatize the reformulation (both exact and approximative) of the models in order to ease the solving part. This is also a promising extension we intend to perform in order to improve significantly the efficiency of this approach. In particular, we intend to rely on [START_REF] Liberti | Reformulation techniques in mathematical programming[END_REF] to implement the reformulations.

Appendix

Drosophila model in AMPL

We present in this appendix how the model presented in Section 3 is written in AMPL. One may see in particular that the formulation is very close to the biological model as it is proposed in [START_REF] Jaeger | Dynamical analysis of regulatory interactions in the gap gene system of drosophila melanogaster[END_REF][START_REF] Nanfack | Efficient parameter estimation for models of pattern formation in early embryogenesis of Drosophila melanogaster[END_REF].

• Sets:

1. set N ι ⊂ N of nuclei; 2. set N γ ⊂ N of genes; 3. set T ⊂ N of time instants.
4. set T data ⊆ T of time instants which we have observed data for.

• Parameters:

1. ∆t ∈ R: interval length for time discretization; 2. for i ∈ N ι , a ∈ N γ let g data a i (t) ∈ R |T data |
be the observed data: for t ∈ T data , g data a i (t) is the observed concentration level of gene a in nucleus i at time t; 3. for a ∈ N γ , v a max is the maximum value for gene a 4. v bcd max is the maximum value for gene bcd; 5. for i ∈ N ι , t ∈ T data , g data bcd i (t) ∈ R is the observed data with respect to gene bcd; 6. λ L , λ U : lower and upper bounds for λ; 9. for all a ∈ N γ , i ∈ N ι , init a i ∈ R is the concentration level of gene a at the first time instant (boundary conditions for time=0).

• Variables: For all gene a ∈ N γ , we define:

1. for all b ∈ N γ such that (b, a) ∈ A, W a b is the weight on (b, a); 2. m a is the influence of gene bcd on gene a;

3. h a is the threshold parameter for the sigmoid regulation function; 4. λ a is the decay coefficient for gene a;

5. R a is the production rate for gene a;

6. D a is the diffusion coefficient for gene a;

7. for all i ∈ N ι , t ∈ T g i a (t) ∈ R represents the estimated data: g a i (t) is the concentration level of gene a in nucleus i at time t;

8. for all i ∈ N ι , t ∈ T , u a i (t) ∈ R is an auxiliary function.

∀a ∈ N γ , i ∈ N ι g a i (0) = init a i ; ∀a ∈ N γ    R L ≤ R a ≤ R U λ L ≤ λ a ≤ λ U D L ≤ D a ≤ D U
3. boundary conditions:

∀v ∈ V x 0 v = C 0 (v).

7 .

 7 R L , R U : lower and upper bounds for R; 8. D L , D U : lower and upper bounds for D;

•

  Objective function:

	min	(g a i (t) -g data	a i (t)) 2 +		(W a b v b max ) 2 +	((m a v bcd max ) 2 + h 2 a ).
	a∈N γ i∈N ι				a∈N γ b∈N γ	a∈N γ
	t∈T data
						(11)
	• Constraints:
	1. the discretized differential equation (2): ∀i ∈ N ι {0, |N ι | -1}, t ∈
	T {0}
	g a i (t)-g a i (t -1) = ∆t	R a 2	(	u a i (t) u a i (t) 2 + 1	+ 1) -λ a g a i (t) + D a (g a i+1 (t) -2g a i (t) + g a i-1 (t)) ;
						(12)
	2. the definition of u i :
	∀i ∈ N ι , t ∈ T u a i (t) =	W a b g b i (t) + m a g bcd i (t) + h a ;	(13)
						b∈N γ
	3. the initial conditions and ranges:

Acknowledgments

This work was supported by EU FP6 FET Project "Morphex". We are also deeply grateful to Françoise Monéger and Jan Traas from the RDP-ENS-Lyon institute for having provided the data on arabidopsis on which we tested our approach.

Arabidopsis model in AMPL

• Sets:

1. set V of genes in the network; 2. set A of arcs in the network; 3. set T of time instants. 4. set U ∈ V of genes involved in the stable sub-network.

• Parameters:

is the initial configuration of the network (vector of boolean values assigned to the genes). 2. α : A → {+1, -1} is the sign of the arc weights; 3. θ L , θ U are the bounds on the threshold values; 4. S : U → {0, 1} is the targeted configuration of the sub-network composed by U .

• Variables:

1. for all i ∈ V , t ∈ T , x t i ∈ {0, 1} is the activation status of gene i at time t; 2. s : T → {0, 1} is a decision variable indicating that the sub-network is stable during at least two successive time steps. 3. y : T → {0, 1} is a decision variable that indicates the first time the sub-network reaches a stable state. 4. θ : V → Z is the threshold function; 5. w : A → N is the arc weight function.

• Objective function:

• Constraints:

1. evolution rule:

2. fixed point conditions: