L O Ï C P A U L E V É

M O R G A N M A G N I N, O L I V I E R R O U X

Static Analysis by Abstract Interpretation of Biological Regulatory Networks Dynamics

The analysis of dynamics of Biological Regulatory Networks (BRNs) requires innovative methods to cope with the state space explosion. This paper settles an original approach for deciding reachability properties based on the Process Hitting, a framework suitable to model dynamical complex systems. In particular, the Process Hitting has been shown of interest to model dynamics of BRNs with discrete values. The Process Hitting describes the way each process is able to act upon (i.e. to "hit") an other one (or itself) in order to "bounce" it as another process further acting. By using complementary abstract interpretations of the succession of actions, we build a very efficient static analysis to over-and under-approximate reachability properties within Process Hittings. Applied to a large BRN of 94 components, our method replies quasi-instantaneously to reachability questions, overtaking the state-of-the-art approaches and showing a very promising scalability.

Introduction

Biological regulatory networks (BRNs) are a common framework to model the concurrent regulations between biological components (RNA, proteins, etc.). These regulations are generally represented as interaction graphs, where nodes are components of the system, and edges state the regulation between them, either positive (activation) or negative (inhibition). To each node is also assigned a numerical value representing the state (e.g. the concentration) of the component of the network, at a given time. This value evolves in response to the various regulations the component is subject to. In 1973, the biologist René Thomas proposed a formalisation of BRNs where the value of components are boolean [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF]. His formalisation uses an interaction graph and René Thomas' parameters (or equivalently, boolean functions between nodes inputs) to describe dynamics of a BRN. A full description of BRNs formalism with discrete values for components can be found in [START_REF] Bernot | Semantics of biological regulatory networks[END_REF].

The derivation of dynamical properties from the interaction graph of BRNs has been the motivation of various mathematical works. Twenty years ago, René Thomas conjectured that the presence of positive circuits within the interaction graph is a necessary condition to achieve systems with multi-stationnarity. The conjecture has been proven in several frameworks, notably in discrete dynamical systems [START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF]. By using more elaborated interaction graph analyses, the maximum number of fixed points within boolean networks can be characterised [START_REF] Aracena | Maximum number of fixed points in regulatory boolean networks[END_REF]. Under strong conditions, particular fixed points (qualified as topological) can be fully extracted from the interaction graph [START_REF] Paulevé | Topological Fixed Points in Boolean Networks[END_REF]: these points are fix in all possible dynamics matching the interaction graph. Finally, the presence of negative circuits in interaction graphs has been proven necessary for sustained oscillations in the dynamics [START_REF] Remy | Graphic requirements for multistability and attractive cycles in a boolean dynamical framework[END_REF][START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF].

To produce more precise analyses of BRNs dynamics, it is required to take into account the boolean functions specified together with the interaction graph (as in [START_REF] Bernot | Gene regulatory networks with multiplexes[END_REF], for instance). The majority of current techniques use standard modelchecking methods [START_REF] Richard | Modern Formal Methods and Applications, chapter Formal Methods for Modeling Biological Regulatory Networks[END_REF] based on the (explicit or symbolic) exploration of the state space of the model. Such methods suffer from the state space explosion, and are intractable on large regulatory networks. We propose here a novel and original method relying on the Process Hitting framework to address this scalability issue.

The Process Hitting [START_REF] Paulevé | Refining Dynamics of Gene Regulatory Networks in a Stochastic π-Calculus Framework[END_REF]) is a recently introduced framework suitable to model BRNs with discrete values. Basically, each discrete component value is modelled as a process; at any time, one and only one process of each component (referred to as sort) is present; this process stands for the current state of the component. A sort changes of process on the hit of at most one other process. Static analyses have already been developed in the Process Hitting framework, notably for obtaining all the fixed points of dynamics of a Process Hitting [START_REF] Paulevé | Refining Dynamics of Gene Regulatory Networks in a Stochastic π-Calculus Framework[END_REF]. Being a particular restriction of Communicating Finite-State Machines [START_REF] Brand | On communicating finite-state machines[END_REF], the Process Hitting can be applied to less specific dynamical complex systems.

The static analysis by abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] aims at providing an efficient analysis of a program without executing it. This is achieved by constructing one or several sound abstractions of the semantics of the program; these abstractions are then interpreted to decide the validity of a given property, resulting in over-or under-approximations of the validity of the property in the concrete program. Hereafter, we refer to this validity as the concretizability of the abstract property.

In this paper, we present a novel static analysis by abstract interpretation of Process Hittings. We address the decision of a successive reachability of processes. Our approach is based on two complementary abstractions of a succession of actions within a Process Hitting. Several refinement operators upon these two abstractions are then defined. These refinements detail an abstraction, with the aim of simplifying the concretizability decision. By using the abstraction refinement operators, over-and under-approximations of the process reachability decision are developed. Their implementations rely on the analysis of an abstract structure, that can be represented as a graph. We show that this abstract structure has always a reasonable size (i.e., polynomial in the total number of processes); and, while its computation can be exponential in the number of processes within a single sort, the approximations are always linear or polynomial in its size.

The scalability of our approach is illustrated by its application to the decision of reachability of gene expression levels within a BRN of 94 components. Our methods responds very fast to the decision, while a well established symbolic model checking technique (SDD, [START_REF] Hamez | Building efficient model checkers using hierarchical set decision diagrams and automatic saturation[END_REF]) regularly fails because of the state space explosion.

A preliminary version of these results have been presented in [START_REF] Paulevé | Abstract Interpretation of Dynamics of Biological Regulatory Networks, The First International Workshop on Static Analysis and Systems Biology[END_REF]. This paper extends them in several ways: the general framework for presenting the analyses has been deeply reworked and unified; the conclusiveness of the overall method has been largely improved; the BRN application case has been extended to 94 components, instead of 40.

This paper is structured as follows. The Process Hitting framework is defined in Sect. 2. Sect. 3 presents complementary abstractions of scenarios (i.e., sequences of actions) and defines the abstraction refinements operators. These refinements operators are then applied to the over-and under-approximation of process reachability in Sect. 4; their implementation and complexity is also detailed. Sect. 5 briefly presents the encoding of BRNs dynamics into Process Hittings and applies the above methods to a large BRN relating 94 components. Finally, Sect. 6 summarises and discusses the contributions of this paper.

Notations: Given a countable set S = {e 1 , . . . , e n }, |S| = n; ℘(S) is the power set. Given a finite sequence of elements A = e 1 :: . . . :: e n , |A| = n is the length of the sequence;

I A = {1, . . . , |A|} is the set of A indexes; A i = e i , ∀i ∈ I A ; ε is the empty sequence; A i..j is the subsequence A i , . . . , A j ; A i..j = ε if j < i. lfp{x 0 } (x → x ) is the least fix point of the function f (x) = x initially applied to x 0 .

The Process Hitting Framework

This section presents the Process Hitting framework [START_REF] Paulevé | Refining Dynamics of Gene Regulatory Networks in a Stochastic π-Calculus Framework[END_REF] on which the methods presented in this paper rely.

The Process Hitting gathers a finite number of concurrent processes grouped into a finite set of sorts. A process belongs to one and only one sort and is noted a i where a is the sort and i the identifier of the process within the sort a. At any time, one and only one process of each sort is present, forming a state of the Process Hitting.

The concurrent interactions between processes are defined by a set of actions. Actions describe the replacement of a process by another of the same sort conditioned by the presence of at most one other process in the current state of the Process Hitting. An action is denoted by

a i → b j b k where a i , b j , b k are processes of sorts a and b. It is required that b j = b k and that a = b ⇒ a i = b j . An action h = a i → b j b k is
read as "a i hits b j to make it bounce to b k ", and a i , b j , b k are called respectively hitter, target and bounce of the action, and can be referred to as hitter(h), target(h), bounce(h), respectively.

Definition 1 (Process Hitting). A Process Hitting is a triple (Σ, L, H):

-Σ = {a, b, . . . } is the finite countable set of sorts, -L = a∈Σ L a is the set of states, with L a = {a 0 . . . a la } the finite and countable set of processes of sort a ∈ Σ and l a a positive integer,

a = b ⇒ a i = b j ∀(a i , b j ) ∈ L a × L b , -H = {a i → b j b k , • • • | (a, b) ∈ Σ 2 , (a i , b j , b k ) ∈ L a ×L b ×L b , b j = b k , a = b ⇒ a i = b j },
is the finite set of actions.

Proc refers to the set of all processes (Proc =

{a i | a ∈ Σ ∧ a i ∈ L a }).
The sort of a process a i is referred to as Σ(a i ) = a and the set of sorts present in an action h ∈ H as Σ(h) = {Σ(hitter(h)), Σ(target(h))}. Given a state s ∈ L, the process of sort a ∈ Σ present in s is denoted by s[a], that is the a-coordinate of the state s. We define the following notations: if

a i ∈ L a , a i ∈ s ⇔ s[a] = a i ; and if ps ∈ ℘(P rocs), ps ⊆ s ⇔ ∀a i ∈ ps, a i ∈ s. An action h = a i → b j b k ∈ H is playable in s ∈ L if and only if s[a] = a i and s[b] = b j .
In such a case, (s • h) stands for the state resulting from the play of the action h in s, that is (s

• h)[b] = b k and ∀c ∈ Σ, c = b, (s • h)[c] = s[c]. For the sake of clarity, ((s • h) • h ), h ∈ H is abbreviated as (s • h • h ).
If A is a sequence of actions, the set of sorts present in A is given by Σ

(A) = n∈I A Σ(A n ).
The first (resp. last) process of sort a appearing in the sequence is referred to as fst a (A) (resp. last a (A)):

fst a (A) =        ∅ if a / ∈ Σ(A), hitter(A m ) if m = min{n ∈ I A | a ∈ Σ(A n )} ∧ Σ(hitter(A m )) = a, target(A m ) else if Σ(target(A m )) = a ;
(1)

last a (A) =            ∅ if a / ∈ Σ(A), bounce(A m ) if m = max{n ∈ I A | a ∈ Σ(A n )} ∧Σ(bounce(A m )) = a, hitter(A m ) else if Σ(hitter(A m )) = a .
(2)

Amongst sequences of actions, the particular sequences composed only of successively playable actions form scenarios (Def. 2). A scenario δ is said to be playable in a state s ∈ L, if and only if δ 1 is playable in s and for all n ∈ I δ , n < |δ|, δ n+1 is playable in the state (s•δ 1 •. . .•δ n ); or equivalently δ is playable in s if and only if support(δ) ⊆ s (Eq. ( 3)). The state resulting from the sequential play of the scenario in s is denoted by s • δ. One can easily show that ∀a i ∈ end(δ), (s

• δ)[a] = a i and ∀b ∈ Σ \ Σ(δ), (s • δ)[b] = s[b];
where end(δ) is defined in Eq. (4).

Definition 2 (Scenario (Sce)). Given a Process Hitting (Σ, L, H), a scenario δ is a sequence of actions in H such that for all n ∈ I δ , a i = hitter(δ n ) (resp. target(δ n )) ⇒ last a (δ 1..n-1 ) ∈ {∅, a i }. The set of all scenarios is denoted by Sce. support(δ) and end(δ) give the first and last processes of each sort, respectively:

support(δ) = {p ∈ Proc | Σ(p) ∈ Σ(δ) ∧ p = fst Σ(p) (δ)} , (3) end
(δ) = {p ∈ Proc | Σ(p) ∈ Σ(δ) ∧ p = last Σ(p) (δ)} . (4)
= {a 0 , a 1 } × {b 0 , b 1 , b 2 } × {c 0 , c 1 } × {d 0 , d 1 , d 2 } and H = {a 0 → c 0 c 1 , a 1 → b 1 b 0 , c 1 → b 0 b 1 , b 1 → a 0 a 1 , b 0 → d 0 d 1 , b 1 → d 1 d 2 , d 1 → b 0 b 2 , c 1 → d 1 d 0 , b 2 → d 0 d 2 }. Playing the action b 1 → a 0 a 1 in the state a 0 , b 1 , c 0 , d 0 results in the state a 1 , b 1 , c 0 , d 0 . δ = a 0 → c 0 c 1 :: b 1 → a 0 a 1 :: a 1 → b 1 b 0 :: b 0 → d 0 d 1 :: d 1 → b 0 b 2 is a scenario playable in the state s = a 0 , b 1 , c 0 , d 0 , and s • δ = a 1 , b 2 , c 1 , d 1 .
Remark 1. The Process Hitting framework can be considered as a class of Communicating Finite-State Machines [START_REF] Brand | On communicating finite-state machines[END_REF] where at most two machines (sorts) share a synchronization label (action) and one and only one machine changes its state (process) at each synchronization (action play).

Abstract Interpretation of Scenarios

After having introduced preliminary definitions (Subsect. 3.1), this section establishes two orthogonal abstractions of scenarios: by objective sequences (Subsect. 3.2) and by bounce sequences (Subsect. 3.3). The former describes a succession of process changes per sort (called objectives), while the latter details the actions actually played to resolve these objectives. While objective sequences can be seen as a sparse representation of a scenario, bounce sequences emphases the necessary actions required to resolve an objective.

Upon these two complementary abstractions, several objective sequence refinements operators are derived in Subsect. 3.4. The aim of these refinements is at providing more precise abstractions on which the concretizability may be easier to decide. Fig. 2 summarizes the possible derivations between a Process Hitting and the different representations of scenarios.

Preliminaries

This subsection introduces the notion of objective and context used by the developed abstractions.

The reachability of the process a j from a process a i is called an objective and is denoted by a i * a j (Def. 3).

Definition 3 (Objective (Obj)). The reachability of process a j from a i is called an objective, noted a i * a j . The set of objectives is denoted by

Obj = {a i * a j | a ∈ Σ ∧ (a i , a j ) ∈ L 2 a }.
Given an objective P ∈ Obj, where P = a i * a j , Σ(P ) = a, target(P ) = a i , bounce(P ) = a j . An objective P is trivial if target(P ) = bounce(P ).

We extend the notion of state by the notion of context (Def. 4). A context references the set of processes per sort that can serve as initial state.

Definition 4 (Context ς (Ctx)). A context ς associates with each sort in Σ a nonempty subset of its processes:

∀a ∈ Σ, ς[a] ⊆ L a ∧ ς[a] = ∅. Ctx refers to the set of contexts. Given a context ς, we note a i ∈ ς ⇔ a i ∈ ς[a]; ps ∈ ℘(Proc), ps ⊆ ς ⇔ ∀a i ∈ ps, a i ∈ ς.
The override of a context ς by a set of processes ps is noted ς ps (Def. 5). For instance,

a 1 , a 2 , b 1 , c 1 {a 3 , b 2 , b 3 } = a 3 , b 2 , b 3 , c 1 .
Definition 5 ( : Ctx×℘(Proc) → Ctx). Given a context ς ∈ Ctx and ps ∈ ℘(Proc), the override of ς by ps is noted ς ps and is defined as

∀a ∈ Σ, (ς ps)[a] = {p ∈ ps | Σ(p) = a} if ∃p ∈ ps, Σ(p) = a, ς[a] otherwise.
A scenario δ ∈ Sce is playable in the context ς if and only if support(δ) ⊆ ς. The play of δ in ς is denoted by ς • δ where ς • δ = ς end(δ).

Abstraction of Scenarios into Objective Sequences

During the execution of a scenario, processes of different sorts bounce one after the other, following the play of the actions. An abstraction of such an execution is a succession of objectives: the process a j is reached (after a certain number of actions) from a i , then the process b j is reached from b i , etc. This forms an objective sequence (Def. 6). The append of an objective to an objective sequence is specified in Def. 7.

Definition 6 (Objective Sequence (OS)). An Objective Sequence is a sequence ω = P 1 :: . . . :: P |ω| , where ∀n ∈ I ω , ω n ∈ Obj and

a i = target(ω n ) ⇒ last a (ω 1..n-1 ) ∈ {∅, a i }.
The set of objective sequences is referred to as OS. The definitions of last a (Eq. ( 2)), fst a (Eq. ( 1)), support (Eq. ( 3)) and end (Eq. ( 4)) are straightforwardly derived to objective sequences by omitting the hitter case.

Definition 7 (⊕ : OS × Obj → OS). The join between an objective sequence ω ∈ OS and an objective a i * a j ∈ Obj is defined as:

ω ⊕ a i * a j = ω :: a i * a j if a / ∈ Σ(ω),
ω :: last a (ω) * a j otherwise.

A scenario can be abstracted by several objective sequences, describing process changes more or less sparsely. For instance, the scenario

a 0 → c 0 c 1 :: b 1 → a 0 a 1 :: a 1 → b 1 b 0 :: b 0 → d 0 d 1 :: d 1 → b 0 b 2 can be abstracted to c 0 * c 1 :: a 0 * a 1 :: b 1 * b 0 :: d 0 * d 1 :: b 0 * b 2 ;
or, more sparsely, to a 0 * a 1 :: b 1 * b 2 ; or to b 1 * b 2 ; etc. (bolded processes are the ones kept in the succeeding abstraction).

The set of scenarios concretizing an objective sequence ω in a context ς is given by γ ς (ω) (Def. 8). We also define the concretization of a set of objective sequences as the union of their concretizations (Def. 9).

Definition 8 (γ

ς : OS → ℘(Sce)). Given ω ∈ OS, γ ς (ω) is the set of scenarios concret- izing ω in the context ς: γ ς (ω) = {δ ∈ Sce |(ω = ε ∧ δ = ε) ∨ ω = ε ∧ support(δ) ⊆ ς ∧ ∃φ : I ω → I δ , (∀n, m ∈ I ω , n < m ⇔ φ(n) ≤ φ(m)) ∧ ∀n ∈ I ω , bounce(ω n ) ∈ ς • δ 1..φ(n) } ,
where ω refers to the objective sequence ω where trivial objectives have been removed.

Definition 9 (γ ς : ℘(OS) → ℘(Sce)). γ ς (Ω) = {δ ∈ γ ς (ω) | ω ∈ Ω} .
It is worth noticing that the concretization of an objective sequence ω does not depend on its support support(ω) as it is imposed by the context. In that way, we use * a i to denote an objective where the target can be any process of sort a present in the context:

γ ς ( * a i ) = γ ς (a j * a i ), ∀a j ∈ ς[a] . (5)
We finally refer to α ς as the reverse operation of concretization of a set of scenarios (Def. 10). This straightforwardly provides Property 1.

Definition 10 (α ς : ℘(Sce) → ℘(OS)). α ς (∆) = {ω ∈ OS | ∃δ ∈ ∆, δ ∈ γ ς (ω)} . Property 1. ∀∆ ∈ ℘(Sce), ∆ ⊆ γ ς (α ς (∆)) . Proof. ∀δ ∈ ∆, there exists ω ∈ OS such that δ ∈ γ ς (ω) (Def. 8). Hence, ω ∈ α ς (∆) (Def. 10), therefore δ ∈ γ ς (α ς (∆)).

Abstraction of Scenarios into Bounce Sequences

Bounce sequences result from a local reasoning on a single sort a. Bouncing from a i to a j (i.e., resolving the objective a i * a j ) may require the play of several actions on processes of sort a, forming a bounce sequence (Def. 11). Remark that bounce sequences are generally not scenarios: e.g., b i → a i a j ::

b j → a j a k is a bounce sequence but not a scenario if b i = b j . Definition 11 (Bounce Sequence (BS)). A bounce sequence ζ is a sequence of actions such that ∀n ∈ I ζ , n < |ζ|, bounce(ζ n ) = target(ζ n+1
). BS denotes the set of bounce sequences. We refer to the set of bounce sequences resolving the objective P as BS(P ):

BS(a i * a j ) = {ζ ∈ BS | target(ζ 1 ) = a i ∧ bounce(ζ |ζ| ) = a j } . Obviously, BS(a i * a i ) = {ε}; and BS(a i * a j ) = ∅ if there is no possibility to reach a j from a i .
In a bounce sequence ζ, target and bounces of all actions share the same sort Σ(ζ). In the scope of this paper, we do not consider bounce sequences that contain cycles between targets and bounces of actions. In that way, the maximum length of a bounce sequence for a sort a is the number of processes of sort a.

The full set of bounce sequences can be computed directly from the set of actions H of the Process Hitting without any enumeration of scenarios. Given an objective a i * a j , the computation of bounce sequences BS(a i * a j ) (Def. 11) works by a depth-first research between actions on the sort a to form a bounce sequence without cycle. Such a computation is exponential in the number of actions on the sort a, and is then efficient when this number is small in front of the total number of actions.

We also consider a sparse representation of a bounce sequence ζ resolving an objective P by considering the set of hitters of its actions that have a different sort than that of P . We denote by BS ∧ (P ) the set of such abstracted bounce sequences (Def. 12).

Definition 12 (BS ∧ : Obj → ℘(Proc)).

BS ∧ (P ) = {ζ ∧ | ζ ∈ BS(P ), ζ ∈ BS(P ), ζ ∧ ζ ∧ } , where ζ ∧ = {hitter(ζ n ) | n ∈ I ζ ∧ Σ(hitter(ζ n )) = Σ(P )}.
It is worth noticing that BS ∧ (P ) can be computed directly from the Process Hitting in the same manner as BS(P ), but yet more efficiently since only minimal sets of hitters are kept, pruning redundant explorations.

The relations of abstractions and concretizations between scenarios and (abstracted) bounce sequences can be derived easily, and are not detailed here.

Looking at the Process Hitting example in Fig. 1,

ζ = a 1 → b 1 b 0 :: d 1 → b 0 b 2 is the only bounce sequence resolving the objective b 1 * b 2 (i.e. BS(b 1 * b 2 ) = {ζ}, and BS ∧ (b 1 * b 2 ) = {{a 1 , d 1 }}).

Objective Sequence Refinements

Before introducing objective sequence refinements operators, we define the relation OS between two objective sequences (Def. 13). Basically, if ω OS ω , we say ω is a more precise abstraction than ω , hence γ ς (ω) ⊆ γ ς (ω ) (Property 2). In such a setting, an objective sequence joined to an objective is a more precise abstraction than the objective alone (Property 3).

Definition 13 ( OS ⊂ OS × OS). ω OS ω if and only if the following properties are satisfied:

-|ω| ≥ |ω |; there exists a mapping φ :

I ω → I ω such that ∀n ∈ I ω , bounce(ω n ) = bounce(ω φ(n) ) and ∀n, m ∈ I ω , n < m ⇔ φ(n) < φ(m). Property 2. ω OS ω =⇒ γ ς (ω) ⊆ γ ς (ω ) .
Property 3. Given ω ∈ OS and P ∈ Obj, ω ⊕ P OS P .

Finally, given an objective P , BS ς (P ) (Def. 14) and BS ∧ ς (P ) (Def. 15) generalise BS(P ) and BS ∧ (P ) to the scope of the context ς, respectively. Finally, Lemma 1 states the correctness of this refinement, ensuring the preservation of the concretization set.

Definition 14 (BS ς : Obj → ℘(BS)). BS ς ( * a j ) = ai∈ς[a] BS(a i * a j ) . Definition 15 (BS ∧ ς : Obj → ℘(℘(Proc))). BS ∧ ς ( * a j ) = ai∈ς[a] BS ∧ (a i * a j ) . 3.
Definition 16 (β : BS → ℘(OS)).

β(ζ) = {ω ∈ OS | |ω| = |ζ| ∧ ∀n ∈ I ζ , bounce(ω n ) = hitter(ζ n )} . Definition 17 (ρ : Obj × ℘(BS) → ℘(OS)). ρ(P, zs) = {ω ⊕ P | ω ∈ β(ζ), ζ ∈ zs} .
Lemma 1. γ ς (P ) = γ ς (ρ(P, BS ς (P ))) .

Proof. (⊇) ∀ω ∈ ρ(P, BS ς (P )), ω OS P , therefore γ ς (P ) ⊇ γ ς (ρ(P, BS ς (P ))); (⊆) By definition of BS ς (P ), ∀δ ∈ γ ς (P ), ∃ω ∈ ρ(P, BS ς (P )), δ ∈ γ ς (ω), thus γ ς (P ) ⊆ γ ς (ρ(P, BS ς (P ))).

3.4.2. Objective refinement by BS ∧ (ρ ∧ ). The refinement of an objective P by BS ∧ is done in a similar way. A set of hitters ps ∈ BS ∧ (P ) is abstracted by β ∧ (ps) into the set of objective sequences describing any ordering of reach of these hitters. The relation between objective sequences in β(ζ) and in β ∧ (ps) is emphased in Property 4. The refinement ρ ∧ (P, BS ∧ ς (P )) is presented in Def. 19 and the preservation of concretizations is stated in Lemma 2.

Definition 18 (β ∧ : ℘(Proc) → ℘(OS)).

β ∧ (ps) = {ω ∈ OS | |ω| = |ps| ∧ ∀p ∈ ps, ∃n ∈ I ω , bounce(ω n ) = p} . Property 4. ∀ζ ∈ BS(P ), ∀ω ∈ β(ζ), ∃ω ∈ β ∧ (ζ ∧ ), ω OS ω . Definition 19 (ρ ∧ : Obj × ℘(℘(Proc)) → ℘(OS)). ρ ∧ (P, pss) = {ω ⊕ P | ω ∈ β ∧ (ps), ps ∈ pss} . Lemma 2. γ ς (P ) = γ ς (ρ ∧ (P, BS ∧ ς (P ))) .
Proof. (⊇) ∀ω ∈ ρ ∧ (P, BS ∧ ς (P )), ω OS P ; (⊆) by Property 4, ∀ω ∈ ρ(P, BS ς (P )), ∃ω ∈ ρ ∧ (P, BS ∧ ς (P )), ω OS ω , thus γ ς (ρ(P, BS ς (P )) ⊆ γ ς (ρ ∧ (P, BS ∧ ς (P )).

3.4.3. Objective sequence refinements ( ρ). Finally, to generalize the refinements defined on objective to objective sequence, we exhibit an objective sequence refinement operator that uses any of the above refinements. We took the operator ρ as an example and define the refinement ρ(ω, BS) (Def. 21). Basically, this refinement chooses any objective ω n of the objective sequence, refines it using ρ, and returns all the interleaving of the obtained refined sequences with the objective sequence ω 1..n-1 (Def. 20); the concretization set is then preserved (Lemma 3).

Definition 20 (interleave

: OS × OS → ℘(OS)). interleave(ω 1 , ω 2 ) = {ω ∈ OS ||ω| = |ω 1 | + |ω 2 | ∧ ∃φ 1 : I ω 1 → I ω , φ 2 : I ω 2 → I ω , (∀n, m ∈ I ω 1 , n < m ⇔ φ 1 (n) < φ 1 (m)) ∧(∀n, m ∈ I ω 2 , n < m ⇔ φ 2 (n) < φ 2 (m)) ∧( n 1 ∈ I ω 1 , n 2 ∈ I ω 2 , φ 1 (n 1 ) = φ 2 (n 2 ))} . Definition 21 ( ρ : OS × ℘(BS) → ℘(OS)). ρ(ω, BS) = { ⊕ ω n..|ω| |n ∈ I ω , ω ⊕ ω n ∈ ρ(ω n , BS ς (ω n )), ∈ interleave(ω 1..n-1 , ω )} . Lemma 3. γ ς (ω) = γ ς ( ρ(ω, BS)) . Proof. (⊇) ∀ω ∈ ρ(ω, BS), ω OS ω; (⊆) by Lemma 1 and Def. 20, δ ∈ γ ς (ω) ⇒ ∃ω ∈ ρ(ω, BS), δ ∈ γ ς (ω ).

Over-and Under-approximations of Process Reachability

We define the Process Reachability problem as deciding if a given objective sequence ω ∈ OS is concretizable for a given Process Hitting in a context ς; i.e. if the set γ ς (ω) is not empty. The process reachability problem can also be formulated in a subclass of the CT L [START_REF] Clarke | Design and synthesis of synchronization skeletons using branching-time temporal logic, Logic of Programs[END_REF], restricted to the following form:

Φ ::= a i | a i ∧ ϕ ϕ ::= EF Φ ,
where a i ∈ Proc is true if the current state contains the process a i and EF stands for the usual exists finally predicate. Given an objective sequence ω, one can encode it into CTL using the following recursive definition of [| • |]:

[|a j * a i :: ε|] = EF a i [|a j * a i :: ω|] = EF (a i ∧ [|ω|]) if ω = ε .
Based upon the refinements operators defined in the previous section, we establish several necessary or sufficient conditions for the concretizability of an objective sequence in a given context ς. These objective sequences can be either given by a user (to check some temporal properties), or extracted from BS (with β, Def. 16) to refine this set of bounce sequences. Indeed, if a bounce sequence is not concretizable in ς, it can be ignored in all analyses in the scope of this context.

These approximations aim at being very fast to compute, overcoming the state space explosion problem inherent to such a kind of dynamics analysis. While inconclusive in some cases, the application section (Sect. 5) shows the very good suitability of our analyses to biological regulatory networks dynamics with a very promising scalability.

The remainder of this section is structured as follows. Subsect. 4.1 presents a first overapproximation based on an un-ordered analysis of objectives required to concretize the given objective sequence; Subsect. 4.2 refines this former approximation by exploiting the sequentiality of objectives to concretize; then Subsect. 4.3 uses order constraints between process occurrences to complete these over-approximations. Finally, Subsect. 4.4 sets up an under-approximation of the process reachability decision.

Un-ordered Over-approximation

Given a context ς and an objective sequence ω, we obtain that ω is concretizable only if each objective ω n , n ∈ I ω , is independently concretizable in the same context (Proposition 1). In this way, one can approximate the concretizability of a bounce sequence by recursively applying Proposition 1 to extract objectives from the given objective sequence, and use the refinement operator ρ ∧ to extend the objective into several objective sequences.

Proposition 1. γ ς (ω) = ∅ =⇒ ∀n ∈ I ω , γ ς (ω n ) = ∅ .
Proof. By Def. 13 and Property 2, ω OS ω i .

Given an objective P , minCont ς (P ) (Def. 22) is the set of re-targeted objectives Q, with target(P ) = target(Q) and bounce(P ) = bounce(Q), which are always derived from a recursive application of ρ ∧ (P, BS ∧ (P )) with Proposition 1 (Lemma 4). The recursive procedure to check necessary conditions for the concretizability of an objective is then summarised by Proposition 2. If P is concretizable, then there is an execution of this procedure without cycle, and all tested objectives have at least one solution (Theorem 1).

Definition 22 (minCont

ς : Obj → ℘(Obj)). minCont ς ( * a j ) = {a k * a j | a k = a j ∧ ∀a i ∈ ς[a], a k ∈ minCont Obj ς (a, a i * a j )} minCont Obj ς : Σ × Obj → ℘(Proc) minCont Obj ς (a, P ) = ∅ if BS ∧ (P ) = ∅, otherwise, minCont Obj ς (a, P ) = {p ∈ Proc | ∀ps ∈ BS ∧ (P ), ∃q ∈ ps, p ∈ minCont P roc ς (a, q)} minCont P roc ς : Σ × Proc → ℘(Proc) minCont P roc ς (a, b i ) =        {b i } if a = b, {p ∈ Proc | ∀b j ∈ ς[b], p ∈ minCont Obj ς (a, b j * b i )} otherwise. Lemma 4. a k * a j ∈ minCont ς ( * a j ) =⇒ γ ς ( * a j ) = γ ς ( * a k :: a k * a j ) .
Proof. By induction on minCont ς , a k occurs in all recursive refinements of * a j by ρ ∧ .

Proposition 2. γ ς (P ) = ∅ =⇒ ∃ps ∈ BS ∧ ς (P ), ∀p ∈ ps, γ ς ( * p) = ∅, and ∀Q ∈ minCont ς (P ), γ ς (Q) = ∅, with BS ∧ ς ( * a j ) = {ps ∈ BS ∧ (a i * a j ) | a i ∈ ς[a]}.
Proof. By Lemma 2, Lemma 4 and Proposition 1.

Theorem 1 (Un-ordered over-approximation). γ ς (ω) = ∅ =⇒ ∀n ∈ I ω , there exists a finite recursive application of Proposition 2 to γ ς (ω n ) = ∅ such that for all tested objective P , BS ∧ ς (P ) = ∅.

Proof. By induction, if γ ς (P ) = ∅ requires γ ς (P ) = ∅, P = P , and if γ ς (P ) = ∅ requires again γ ς (P ) = ∅, then γ ς (P ) = ∅; therefore, by Proposition 1, there exists a finite recursion of Proposition 2 application on γ ς (ω n ) = ∅, ∀n ∈ I ω . Finally, Proposition 2 implies the non-emptyness of BS ∧ ς (P ) for each tested objective P .

Implementation. Given a context ς and an objective sequence ω, we define an abstract structure A ω ς (Def. 23) which mimics the relations between objectives during the execution of Proposition 2 (Lemma 5). A ω ς gathers three relations (Req ω ς , Sol ω ς , Cont ω ς ), respectively the requirements, solutions, and minimal continuity (or re-targeting).

Definition 23 (A ω ς ). Given a context ς and an objective sequence ω, we define the abstract structure A ω ς = (Req ω ς , Sol ω ς , Cont ω ς ), where Req ω ς , Sol ω ς and Cont ω ς are defined as follows:

Req ω ς ={(a i , a j * a i ) ∈ Proc × Obj | a j ∈ ς[a] ∧ (∃(P, ps) ∈ Sol ω ς , a i ∈ ps ∨ ∃n ∈ I ω , bounce(ω) = a i )} Sol ω ς ={(P, ps) ∈ Obj × ℘(Proc) | ∃(a i , P ) ∈ Req ω ς ∧ ps ∈ BS ∧ (P )} Cont ω ς ={(P, Q) ∈ Obj × Obj | ∃(P, ps) ∈ Sol ω ς ∧ Q ∈ minCont ς (p)} .
Lemma 5. Given an objective P referenced in A ω ς ,

ω ⊕ P ∈ ρ ∧ (P, BS ∧ (P ))

⇐⇒ (P, {bounce(ω n ) | n ∈ I ω }) ∈ Sol ω ς ∧ ∀n ∈ I ω , a j = bounce(ω n ), ∀a i ∈ ς[a], (a j , a i * a j ) ∈ Req ω ς ; -Q ∈ minCont ς (P ) ⇐⇒ (P, Q) ∈ Cont ω ς . Proof. By construction of A ω ς . A ω
ς has a graph structure, with cycles, potentially. We remark that, as |Obj| = a∈Σ |L a | 2 , the size of Req ω ς and Cont ω ς sets are polynomial in the number of processes in the Process Hitting. The size of Sol ω ς also depends on the cardinality of BS ∧ which follows the maximum number of combinations of |L a | different processes, far below an exponential growth.

Finally, Algorithm 1 details the computation of the Theorem 1 decision.

Algorithm 1 (Un-ordered over-approximation). Given a context ς, an objective sequence ω ∈ OS and the abstract structure

A ω ς : 1 Initialise Θ = {P ∈ Obj | (P, ∅) ∈ Sol ω ς }. 2 Repeat until fix-point: (a) Υ = {p ∈ Proc | ∃P ∈ Θ, (p, P ) ∈ Req ω ς }; (b) Θ = {P ∈ Obj | ∃(P, ps) ∈ Sol ω ς , ps ⊆ Υ ∧ ∀(P, Q) ∈ Cont ω ς , Q ∈ Θ}. 3 γ ς (ω) = ∅ =⇒ ∀n ∈ I ω , ∃P ∈ Θ, target(P ) ∈ ς ∧ bounce(P ) = bounce(ω n ).

Complexity. The computation of A ω

ς is done by iteratively adding the required processes and objectives. The steps of Algorithm 1 are done polynomially in the size of A ω ς . Putting aside the BS ∧ computation complexity, the proposed over-approximation can then be achieved by a number of operations polynomial in the size of the abstract structure.

Examples. Fig. 3 represents graphically the abstract structure that has been extracted from the Process Hitting example in Fig. 1 for a particular objective sequence and context. In this case, Theorem 1 is satisfied. Fig. 4 applies Theorem 1 to the Process Hitting example in Fig. 1 for the decision of the objective d 1 * d 2 concretizability. In this case, the necessary condition is not satisfied.

Ordered Over-approximation

This subsection exploits the ordering of objectives in an objective sequence to increase the conclusiveness of the over-approximation for its concretizability.

Given an objective P which is not trivial in the given context ς, we denote by ends(P ) the set of processes a scenario concretizing P may lead to (Def. 24, Proposition 3). This set is computed from BS(P ) by taking the hitter and bounce of the last action in each bounce sequence.

Definition 24 (ends ς : Obj → ℘(℘(Proc))). Proof. As a i / ∈ ς[a], there exists n ∈ I δ such that bounce(δ n ) = a i ; hence δ 1..n ∈ γ ς ( * a i ) and end(δ n ) ⊆ end(δ 1..n ). By Def. 24, end(δ n ) ∈ ends ς ( * a i ).

ends ς ( * a i ) = {end(h) |∃a j ∈ ς[a], ∃ζ ∈ BS(a j * a i ), ζ = ε ∧ h = ζ |ζ| } . Proposition 3. γ ς ( * a i ) = ∅ ∧ a i / ∈ ς[a] =⇒ ∃δ ∈ γ ς ( * a i ), ∃eps ∈ ends ς ( * a i ) such that eps ⊆ end(δ) .
Given an abstract structure A ω ς ∈ A, procs(A ω ς ) is the set of processes referenced in A ω ς (Def. 25) from which derives Proposition 4. We then define maxprocs ς (ω) (Def. 26) as the set of processes present in the abstract structure having its context saturated (i.e., such that ς procs(A ω ς ) = ς). Note a particular optimisation when ω = P , in which case the process bounce(P ) can be ignored by the context.

Definition 25 (procs

: A → ℘(Proc)). procs((Sol ω ς , Req ω ς , Cont ω ς )) = {p ∈ Proc |∃(P, ps) ∈ Sol ω ς , p ∈ ps ∨ p = target(P ) ∨ (P = ω ⇒ p = bounce(P ))} .
Proposition 4. For each objective Q = ω tested by Proposition 2 during the checking of Theorem 1 on γ ς (ω), bounce(Q) ∈ procs(A ω ς ) and ς[Σ(Q)] ⊆ procs(A ω ς ). Proof. By Lemma 5.

Definition 26 (maxprocs ς : OS → ℘(Proc)).

maxprocs ς (ω) = procs( A ω ς ) where A ω ς = lfp{A ω ς } A ω ς → A ω ς procs(A ω ς )
.

By intersecting maxprocs ς (ω) and ends ς (ω 1 ), we obtain a context into which ω 2..|ω| concretizability should satisfy Theorem 1, if ω is concretizable in ς. This is stated by Theorem 2, which gives a straightforward refinement of Theorem 1 by taking the sequentiality of objectives in ω into account.

Theorem 2 (Ordered over-approximation). Given a context ς and ω = P :: ω ∈ OS such that bounce(P ) / ∈ ς, γ ς (P :: ω ) = ∅ =⇒ Theorem 1 is satisfied with γ maxς (ω ) = ∅, where maxς = ς maxprocs ς (ω) eps and eps ∈ ends ς (P ).

Proof. δ ∈ γ ς (P :: ω ) ⇒ ∃n ∈ I δ such that bounce(δ n ) = bounce(P ), eps ⊆ end(δ 1..n ), and δ n+1..|δ| ∈ γ ς (ω ) with ς = support(δ). Therefore, Theorem 1 is satisfied with γ ς (ω ) = ∅. Let Q be an objective tested by Proposition 2 when checking Theorem 1 satisfaction with γ maxς (ω ) = ∅. By Proposition 2, target(Q) ∈ maxς. If target(Q) ∈ eps then target(Q) ∈ ς , thus by hypothesis, Proposition 2 is satisfied. In the case of target(Q) / ∈ eps, by Def. 26 and Proposition 4, bounce(Q) ∈ maxprocs ς (ω), thus Proposition 2 holds.

By defining maxCtx(ς, ω, n) as the maximum context after the resolution of ω 1..n (Def. 27), the above theorem can be straightforwardly extended to Corollary 1.

Definition 27 (maxCtx : Ctx × OS × N → Ctx). (we assume n ∈ {0} ∪ I n ): maxCtx(ς, ω, n) =        ς if n = 0, ς maxprocs ς (ω) if bounce(ω n ) ∈ maxCtx(ς, ω, n -1)
ς maxprocs ς (ω) eps otherwise, where eps ∈ ends ς (ω n ) Theorem 2 concludes this objective sequence is not concretizable.

Corollary 1. γ ς (ω) = ∅ =⇒ ∀n ∈ I ω , γ maxς (ω n ) = ∅, with maxς = maxCtx(ς, ω, n -1).
Implementation. The computation of maxprocs ς (ω) requires at most |Proc| iterations, giving a number of steps polynomial in the number of processes. Regarding the computation of ends ς ( * a i ), it can either be derived from prior BS( * a i ) computations; or, if BS are too costly to compute, it can be approximated by either {{bounce(P )}}, or by

{end(h) | h ∈ H ∧ bounce(h) = a i ∧ ∃a j ∈ ς[a], ∃ps ∈ BS ∧ (a j * a i ), target(h) ∈ ps}.
Example. Given the Process Hitting defined in Fig. 5, with ς = a 1 , b 0 , Theorem 1 is inconclusive on the concretizability of ω = a 1 * a 0 :: b 0 * b 1 . By applying Theorem 2, it appears that b 0 * b 1 is not concretizable in maxς = maxCtx(ς, ω, 1) = a 0 , b 0 (in such a case, the A b0 * b1 maxς forms a unique cycle between b 1 and a 1 ).

Over-approximation using Process Occurrences Order Constraints

An objective a i * a j having no solution (BS ∧ (a i * a j ) = ∅) informs that the process a j never occurs after a i . This order constraints between process occurrences is referred to as a j a i (Def. 28, Property 5).

Definition 28 ( ). The binary relation ⊂ Proc×Proc is a partial pre-order such that a j a i (i.e., (a j , a i ) ∈ ) if and only if there exists no scenario where a j occurs after a i :

a j a i ⇐⇒ δ ∈ Sce such that ∃n, m ∈ I δ , n ≤ m, target(δ n ) = a i ∧ bounce(δ m ) = a j . Property 5 (Order constraint uncovering). BS(a i * a j ) = ∅ =⇒ a j a i .
Knowing some order constraints on process occurrences, we want to check if some sequence of objectives doesn't contradict such constraints. This can be achieved by computing the processes that always occur when resolving an objective: given an objective sequence ω, if a process a i is required by ω n and a process a j by ω m , n < m, then the constraint a j a i should not exist. This is illustrated by Fig. 6.

In a similar fashion to minCont ς (Def. 22), minProc ς (P ) refers to the set of processes of any sort that occur in all refinements of P in the context ς (Def. 29, Lemma 6). By . Illustration of the method developed in Subsect. 4.3 to over-approximate the concretizability of an objective sequence ω in the state ς: for each objective, the minimal set of occurring processes (represented as squares) are computed using minProc (Def. 30); the sequence is not concretizable as soon as two processes (in black) occurring in distinct objectives resolution contradict the process occurrences order .

using the previously defined maxCtx (Def. 27), we define minProc(ς, ω, n) (Def. 30) the set of processes occurring in the resolution of ω n after having resolved ω 1..n-1 . From this definition, Theorem 3 states the over-approximation illustrated in Fig. 6.

Definition 29 (minProc ς : Obj → ℘(Proc)). Given a context ς, minProc ς ( * a i ) = {p ∈ Proc | ∀a j ∈ ς[a], BS(a j * a i ) = ∅ ⇒ p ∈ minProc Obj ς (a j * a i )} minProc Obj ς : Obj → ℘(Proc) minProc Obj ς (a j * a i ) = {a i } ∪ {p ∈ Proc | ∀ps ∈ BS ∧ (a j * a i ), ∃q ∈ ps, p ∈ minProc ς ( * q) ∨ ∃a k * a i ∈ minCont Obj ς (a, a j * a i ), p ∈ minProc Obj ς (a i * a k ) ∪ minProc Obj ς (a k * a i )))} .
Lemma 6. ∀δ ∈ γ ς (P ), ∀p ∈ minProc ς (P ), p ∈ δ .

Proof. By induction on minProc ς , p occurs in all recursive refinements of P by ρ ∧ .

Definition 30 (minProc

: Ctx × OS × N → ℘(Proc)). (we assume n ∈ {0} ∪ I ω ): minProc(ς, ω, n) =        {a i ∈ ς} if n = 0 minProc maxς (ω n ) otherwise,
with maxς = maxCtx(ς, w, n -1) .

Theorem 3 (Ordered over-approximation refined by ).

γ ς (ω) = ∅ =⇒ n, m ∈ {0} ∪ I ω , n < m, ∃p ∈ minProc(ς, ω, n), ∃q ∈ minProc(ς, ω, m), q p.
Proof. By Lemma 6, Corollary 1 and Def. 30.

Implementation. The implementation is very similar to the one presented in Subsect. 4.2. The uncovering of is done linearly in the size of the saturated abstract structure A ω ς . Example. Let us define the Process Hitting as in Fig. 7, and its saturated abstract structure A ω ς , with ς = a 1 , b 1 , z 1 and ω = z 1 * z 2 :: z 2 * z 1 :: z 1 * z 2 , for which the concretizability has to be decided. The evaluation of minCont(ς, ω, n) and maxCtx(ς, ω, n) give the following:

n = 0 -ς n = 1 -z 1 * z 2 n = 2 -z 2 * z 1 n = 3 -z 1 * z 2 minProc(ς, ω, n) a 1 , b 1 , z 1 a 0 , a 1 , b 0 , z 0 , z 2 b 1 , z 1 a 0 , a 1 , b 0 , z 0 , z 2 maxCont(ς, ω, n) a 1 , b 1 , z 1 a 0 , a 1 , b 0 , z 2 a 0 , a 1 , b 1 , z 1 -
As a 1 a 0 , ω is not concretizable in ς.

Under-approximation

The under-approximation procedure presented in this subsection takes advantage of a variant of the abstract structure used in the above over-approximations. If certain conditions on this abstract structure are verified, then it is shown that a scenario concretizing the given objective sequence exists. The proposed construction of the scenario is made by a so-called top-down resolution: given an objective sequence ω, we first build the scenario concretizing ω 1 by preempting the resolution of the objective sequence ω 2..|ω| (and hence, ignoring any objective interleaving that may be required to concretize ω).

As stated by the refinement operators in Subsect. 3.4, the concretization of ω 1 involves the concretization of a refinement of ω 1 , resulting in a recursive procedure of scenario construction.

We first define an alternative definition of the set of scenarios concretizing an objective sequence ω in a context ς: ς (ω) is empty unless, for each state s ∈ L included in ς, there exists a scenario δ ∈ γ ς (ω) such that δ is playable in s; in that case, ς (ω) = γ ς (ω) (Def. 31). Property 6 is directly derived from this definition and the extension of ς to a set of objective sequence is given in Def. 32.

Definition 31 ( ς : OS → ℘(Sce)).

ς (ω) = γ ς (ω) if ∀s ∈ L, s ⊆ ς, ∃δ ∈ γ ς (ω), support(δ) ⊆ s ∅ otherwise. Property 6. ς ⊆ ς ∧ ς (ω) = ∅ =⇒ ς (ω) = ∅ . Definition 32 ( ς : ℘(OS) → ℘(Sce)). ς (Ω) = {δ ∈ ς (ω) | ω ∈ Ω}
Given an objective P , maxCont ς (Σ(P ), P ) (Def. 33) is the set of processes of sort Σ(P ) that may be encountered during the resolution of P . We then define the saturated abstract structure B ω ς = ( Req ω ς , Sol ω ς , Cont ω ς ) (Def. 34) similarly to A ω ς (Def. 23), except that Sol ω ς can arbitrarily selects bounce sequences to resolve an objective, and that Cont ω ς reflects maxCont ς instead of minCont ς . It appears that if B ω ς contains no cycle, and if all referenced objectives have at least one solution, then a top-down resolution of any referenced objective succeeds in every state of ς. This is stated by Theorem 4, which provides sufficient conditions for the concretization of an objective sequence in a given context.

Definition 33 (maxCont

ς : Σ × Obj → ℘(Proc)). maxCont ς (a, P ) = {p ∈ Proc |∃ps ∈ BS ∧ (P ), ∃b i ∈ ps, b = a ∧ p = b i ∨ b = a ∧ p ∈ maxCont ς (a, b j * b i ) ∧ b j ∈ ς[b])} . Definition 34 ( B ω ς ). The abstract structure B ω ς = ( Req ω ς , Sol ω ς , Cont ω ς ) is defined as B ω ς = lfp{B ω ς } B ω ς → B ω ς procs(B ω ς ) , with B ω ς = (Req ω ς , Sol ω ς , Cont ω ς ): Req ω ς ={(a i , a j * a i ) ∈ Proc × Obj | a j ∈ ς[a] ∧ (∃(P, ps) ∈ Sol ω ς , a i ∈ ps ∨ ∃n ∈ I ω , bounce(ω) = a i } Sol ω ς ⊆{(P, ps) ∈ Obj × ℘(Proc) | ∃(a i , P ) ∈ Req ω ς ∧ ps ∈ BS ∧ (P )} Cont ω ς ={(P, q * bounce(P )) ∈ Obj × Obj | ∃(P, ps) ∈ Sol ω ς ∧ q ∈ maxCont ς (Σ(P ), P )} .
Theorem 4 (Under-approximation). If the graph B ω ς has no cycle and all referenced objectives have at least one solution, then ς (ω) = ∅. we assume all objectives children of P concretizables (no cycles). If ∃Q ∈ Cont ω ς , then by hypothesis, s (target(P ) * target(Q) :: Q) = ∅, thus s (P ) = ∅. Otherwise, by Def. 33, concretizations of children of P do not require any process of sort Σ(P ). Also, there exists ζ ∈ BS(P ) such that (P,

ζ ∧ ) ∈ Sol ω ς . By hypothesis, ∀n ∈ I ζ , ∃δ n ∈ s n-1 ( * hitter(ζ n )) with either s n-1 = s if n = 1, or s n-1 = s • δ 1 • . . . • δ n-1 ;
and Σ(P ) / ∈ Σ(δ n ) (by Def. 33). Therefore, δ = δ 1 :: ζ 1 :: . . . δ n :: ζ n ∈ s (P ) and end(δ) ⊆ maxς.

Finally, as maxς (ω) = ∅, ς (ω) = ∅ (Property 6).

From the proof of the above theorem, we define endProc( B ω ς , P ) (Def. 35) as the maximum set of processes a scenario built by Theorem 4 may end to (Corollary 2). This allows a straightforward extension of Theorem 4 to take the sequentiality of objectives into account (Corollary 3).

Definition 35 (endProc

: B × Obj → ℘(Proc)). endProc( B ω ς , P ) = {p ∈ Proc | Σ(p) = Σ(P ) ⇒ p = bounce(P ) ∧ (∃(P, Q) ∈ Cont ω ς , p ∈ endProc( B ω ς , target(P ) * target(Q)) ∨ p ∈ endProc( B ω ς , Q) ∨ ∃(P, ps) ∈ Sol ω ς , ∃b i ∈ ps, ∃b j ∈ ς[b], p ∈ endProc( B ω ς , b j * b i ))}
Corollary 2. Given P ∈ Obj, if Theorem 4 is satisfied with ς (P ) = ∅, then ∀s ∈ L, s ⊂ ς, ∃δ ∈ ς (P ) such that end(δ) ⊆ endProc(B P ς , P ).

Corollary 3. Given P ∈ Obj, and ω ∈ OS, if Theorem 4 is satisfied with ς (P ) = ∅, and if Theorem 4 is satisfied with ς (ω) = ∅, where ς = ς endProc(B P ς , P ), then Theorem 4 is satisfied with ς (P ⊕ ω) = ∅.

Implementation. The computation of the saturated abstract structure B ω ς works by progressive addition of relations between objectives and processes, found by several traversing of the abstract structure, giving a maximal complexity polynomial in the size of the abstract structure. Checking for Theorem 4 conditions is done linearly in the size of the obtained abstract structure. It is worth noticing that arbitrarily selecting solutions for objectives in B ω ς prevents spurious saturations and may increase the satisfaction of the above theorem, but potentially increases the complexity of checking (as several combinations of solutions can be tested). Examples. Fig. 8 shows two examples of the application of Theorem 4 on the Process Hitting example in Fig. 1.

Discussion. Corollary 3 suggests that testing a refined objective sequence may reveal more conclusive than testing the original given objective sequence. A future work may use of graph analysis of B ω ς to determine which refined objective sequences are good candidates to satisfy Theorem 4, and then increase the conclusiveness of the method.

Application to Biological Regulatory Networks

From Biological Regulatory Networks to Process Hittings

We first sketch the modelling of a discrete BRN in the Process Hitting framework. Basically, to each component corresponds a sort, and to each state of components corresponds a process. If a component a at state i activates a component b at state j, an action a i → b j b k is added, where b k is the state of b after activation. The inhibition is modelled similarly. The realisation of boolean functions between nodes are modelled using a dedicated sort, and is illustrated in Fig. 9. The full formalisation of this translation can be found in [START_REF] Paulevé | Refining Dynamics of Gene Regulatory Networks in a Stochastic π-Calculus Framework[END_REF].

T-Cell Receptor Signalling Pathway (94 components)

Presented in (Saez-Rodriguez, Simeoni, Lindquist, Hemenway, Bommhardt, Arndt, Haus, Weismantel, Gilles, [START_REF] Saez-Rodriguez | A logical model provides insights into t cell receptor signaling[END_REF], this biological system models the T-Cell Receptor (TCR) signalling pathway, the behaviours of which reveal an activation of transcription factors controlling the cell's fate, e.g. whether it proliferates or not. This model is an extension of a former presented BRN relating 40 components [START_REF] Klamt | A methodology for the structural and functional analysis of signaling and regulatory networks[END_REF]. The Process Hitting model † of this system is composed of 1124 actions between 448 processes splitted in 133 sorts (the largest sort has 16 processes). The total number of states of this model is 2 194 (≈ 2 • 10 58 ).

Independant reachability decisions have been experimented from all possible inputs combinations (components CD45, CD8, TCRlig) to each output (components SRE, AP1, CRE, NFkB, NFAT, Cyc1, p21c p27k, FKHR, BclXL). All result in conclusive decisions. It is worth noticing that only approximations defined by Theorem 1 (Subsect. 4.1) and Theorem 4 (Subsect. 4.4) have been exploited, showing quite simple dynamics for the independant reachability of components.

Computation times are around the hundredth of a second on a 3GHz processor with 2GB of RAM. To give a comparison, we did the same experiments with a standard symbolic model-checking method using state-space compression based on Hierarchical Set Decision Diagrams (SDDs) [START_REF] Hamez | Building efficient model checkers using hierarchical set decision diagrams and automatic saturation[END_REF]: the libddd framework (LIP6/Move, 2007), known for its good performances. For the majority of reachability decisions, the program runs out of memory, for others, computation times range from some seconds to hours. This shows the remarkable efficiency of our method, based on abstract interpretation.

Discussion

The Process Hitting is a recently proposed framework suitable for modelling dynamics of BRNs with discrete values. In Process Hitting, components are represented as sorts, and their levels as processes; at any time, one and only one process of each sort is present. The successive states of a component within the system are enclosed in the so-called sort. The replacement of a process by another of the same sort (i.e. level change of a component), is conditioned by the presence of at most one other process, of any sort.

Thanks to the particular structure of Process Hittings models, a powerful static analysis by abstract interpretation has been developed to decide the successive reachability of processes, hence of component levels in the scope of BRNs modelling. The computation is done by over-and under-approximation of the decision, and may reveal to be inconclusive. It exploits two complementary abstractions of scenarios in Process Hittings (by objective and bounce sequences). Several refinements of an objective sequence are then defined to detail the required steps necessary to its concretizability. These refinements are exploited to derive necessary or sufficient conditions for the process reachability satisfaction. Further work may improve the conclusiveness of our method, as it is discussed in Subsect. 4.4. Also, the link between the conclusiveness of the developed method and the structure of the interaction graph of the BRN could be studied.

The implementation of the presented approximations use of an abstract structure that is ensured to have a size nearly polynomial in the total number of processes. On the one hand, the computation of the refinements can be exponential in the number of processes within a single sort; on the other hand the computations of the decisions are polynomial in the size of the abstract structure. Hence, it is expected efficient analyses when the number of processes within a sort is limited, while a large amount of sorts should be handled.

This new and original approach has been applied to the analysis of a large BRN relating 94 components. Response times are really fast (around the hundredth of a second on a desktop computer), showing a promising scalability of our method. It has been compared to a standard symbolic model checking techniques which regularly fails to analyse the model, because of the state space explosion. To our knowledge, this is the first successful application of model checking to BRNs of such a size.

Related work. There has been recent work on the fast computation the set of reachable components within Kappa models [START_REF] Danos | Abstract interpretation of cellular signalling networks[END_REF] (a rule based language). In the general case, this set is over-approximated and such an analysis does not permit to decide the successive reachability of processes, as it is done in this paper. The work presented in [START_REF] Alimonti | Linear time analysis of properties of conflict-free and general petri nets[END_REF], on so-called T-Paths within Petri Nets, establishes structural properties within Petri Nets to derive either necessary or sufficient conditions for place marking reachability. This follows a similar approach to ours. Finally, [START_REF] Pilegaard | Static analysis of a model of the LDL degradation pathway[END_REF] applies static analysis techniques of bioambient models to study the behaviour of biological systems.

Future work investigate other applications of the presented refinements of the abstraction interpretations of Process Hitting. In particular, as they extract the causality between process changes, they point up processes required for a certain process reachability. This kind of analysis may lead to a control of the studied system by acting upon these key processes: e.g. knockdowning a key gene to prevent a cascade of gene activation (gene therapy).

Another research direction is the incorporation of quantitative aspects within the presented decision of process reachability, such as the probability of reaching a given process in a given time interval.

Fig. 1

 1 Fig. 1 represents a Process Hitting (Σ, L, H) where Σ = {a, b, c, d}, L= {a 0 , a 1 } × {b 0 , b 1 , b 2 } × {c 0 , c 1 } × {d 0 , d 1 , d 2 } and H = {a 0 → c 0 c 1 , a 1 → b 1 b 0 , c 1 → b 0 b 1 , b 1 → a 0 a 1 , b 0 → d 0 d 1 , b 1 → d 1 d 2 , d 1 → b 0 b 2 , c 1 → d 1 d 0 , b 2 → d 0 d 2 }.Playing the action b 1 → a 0 a 1 in the state a 0 , b 1 , c 0 , d 0 results in the state a 1 , b 1 , c 0 , d 0 . δ = a 0 → c 0 c 1 :: b 1 → a 0 a 1 :: a 1 → b 1 b 0 :: b 0 → d 0 d 1 :: d 1 → b 0 b 2 is a scenario playable in the state s = a 0 , b 1 , c 0 , d 0 , and s • δ = a 1 , b 2 , c 1 , d 1 .

Figure 1 Figure 2 .

 12 Figure 1. A Process Hitting example. Sorts are represented by labeled boxes, and processes by circles (ticks are the identifiers of the processes within the sort, for instance, a 0 is the process ticked 0 in the box a). An action (for instance a 0 → c 0 c 1 ) is represented by a pair of directed arcs, having the hit part (a 0 to c 0 ) in plain line and the bounce part (c 0 to c 1 ) in dotted line. The reachability of the process d 2 (double circled) is studied in next sections. The current state is represented by the grayed processes: a 0 , b 1 , c 0 , d 0 .

Figure 3 .Figure 4 .

 34 Figure 3. Graphical representation of the abstract structure A ω ς extracted from the Process Hitting in Fig. 1, with ω = d 0 * d 2 and ς = a 1 , b 1 , b 2 , d 0 .

Figure 5 .

 5 Figure 5. (right) Saturated abstract structure A ω ς of the (left) Process Hitting with ς = a 1 , b 0 and ω = a 1 * a 0 :: b 0 * b 1 . Theorem 2 concludes this objective sequence is not concretizable.

  Figure6. Illustration of the method developed in Subsect. 4.3 to over-approximate the concretizability of an objective sequence ω in the state ς: for each objective, the minimal set of occurring processes (represented as squares) are computed using minProc (Def. 30); the sequence is not concretizable as soon as two processes (in black) occurring in distinct objectives resolution contradict the process occurrences order .

Figure 7 .

 7 Figure 7. (bottom) Saturated abstract structure A ω ς of the (top) Process Hitting with ς = a 1 , b 1 , z 1 and ω = z 1* z 2 :: z 2 * z 1 :: z 1 * z 2 . Theorem 3 concludes this objective sequence is not concretizable.

Proof.

  We denote by maxς = ς procs( B ω ς ) the context handled by B ω ς . By induction on the acyclic graph B ω ς , we prove that ∀s ∈ L, s ⊆ maxς, for each objective P referenced in B ω ς such that targetP ∈ s, ∃δ ∈ s (ω) and end(δ) ⊆ maxς. -(P, ∅) ∈ Sol ω ς ⇒ either target(P ) = bounce(P ) (thus δ = ε), or ∀ζ ∈ BS(P ), ζ ∈ Sce ∧ Σ(ζ) = {Σ(P )}, thus δ = ζ.

Figure 8 .

 8 Figure 8. Saturated abstract structure B ω ς from the Process Hitting in Fig. 1 with ω = d 0 * d 2 and ς = a 1 , b 1 , c 1 , d 0 (top), and ς = a 0 , b 1 , c 0 , d 0 (bottom). By Theorem 4, ω is concretizable in a 1 , b 1 , c 1 , d 0 . At this stage, our procedure is inconclusive for the concretizability of ω in a 0 , b 1 , c 0 , d 0 .

Figure 9 .

 9 Figure 9. Examples of Process Hittings (right) from BRNs, having interaction graphs at left. (1) simple inhibition of c by a. (2) boolean function between a (inhibitor) and b (activator) on c ab reflects the state of sorts a and b. In this case, ab 3 reflects the state a 0 , b 1 (and, ab 0 the state a 1 , b 0 , ab 1 the state a 1 , b 1 , ab 2 the state a 0 , b 0 ).

† Model and implementation available at http://processhitting.wordpress.com


J.-M Friedt

T Rétornaz

S Alzuaga

T Baron

G Martin

T Laroche

S Ballandras

M Griselin

J.-P Simonnet

Surface Acoustic Wave Devices as Passive Buried Sensors

Keywords: numbers: 07.79-v GPR, SAW, buried sensors, passive wireless device

Surface Acoustic Wave (SAW) devices are currently used as passive remote-controlled sensors for measuring various physical quantities through a wireless link. Amongst the two main classes of designs -resonator and delay line -the former has the advantage of providing narrow-band spectrum informations and hence appears compatible with an interrogation strategy complying with ISM (Industry-Scientific-Medical) regulations in Radio-Frequency (RF) bands centered around 434, 866 or 915 MHz. Delay-line based sensors require larger bandwidths as they consists of a few interdigitated electrodes excited by short RF pulses with large instantaneous energy and short response delays, but is compatible with existing equipment such as Ground Penetrating RADAR (GPR). We here demonstrate the measurement of temperature using the two configurations, particularly for long-term monitoring using sensors buried in soil. Although we have demonstrated long term stability and robustness of packaged resonators and signal to noise ratio compatible with the expected application, the interrogation range (max 80 cm) is insufficient for most geology or geophysical purposes. We then focus on the use of delay lines, as the corresponding interrogation method is similar to the one used by GPR which allows for RF penetration distances ranging from a few meters to tens of meters and which operates in the lower RF range, depending on soil water content, permittivity and conductivity. Assuming propagation losses in a pure dielectric medium with negligible conductivity (snow or ice), an interrogation distance of about 40 m is predicted, which overcomes the observed limits met when using interrogation methods specifically developed for wireless SAW sensors, and could partly comply with the above-mentioned applications. Although quite optimistic, this estimate is consistent with the signal to noise ratio observed during an experimental demonstration of the interrogation of a delay line buried at a depth of 5 m in snow.

I. INTRODUCTION

Within the framework of wireless sensors, surface acoustic wave (SAW) piezoelectric devices provide unique performances in terms of robustness and autonomy compared to active devices (better temperature stability compared to Complementary Metal-Oxide-Semiconductor -CMOS -devices, no need for on-board power supply), and larger interrogation distance than Radio-Frequency IDentification -RFID -passive tags. The use of piezoelectric delay lines and resonators for monitoring physical quantities such as temperature, strain, torque and pressure have already been demonstrated [START_REF] Bao | IEEE Ultrasonics Symposium[END_REF][START_REF] Reindl | IEEE Ultrasonics Symposium[END_REF][START_REF] Buff | IEEE Ultrasonics Symposium[END_REF][START_REF] Pohl | [END_REF][5] . As opposed to CMOS based devices, acoustic sensors relying on a piezoelectric substrate do not exhibit an incoming energy level threshold to work properly.

The operating principle of these devices is based on a first conversion of an incoming electromagnetic wave to an acoustic propagating wave. As the latter is sensitive to its environment (and hence can be optimized to sense a specific physical parameter) its principal characteristics, and mainly the phase velocity, are modulated according to the conditions the sensor is submitted to. Finally, the acoustic energy stored in the sensor is converted back to an electrical signal by direct piezoelectric effect and dissipated as an electromagnetic radiation through the antenna. This signal is captured and analyzed to evaluate the measured physical parameter.

The resulting interrogation range and compliance with existing RADAR systems (Ground Penetrating Radar -GPR [START_REF] Daniels | Ground Penetrating Radar[END_REF] ) complements the identification capability already familiar to passive IFF (Indentify Friend/Foe) systems used since the second world war, with measurement capabilities [START_REF] Rieback | [END_REF] .

The short interrogation delay (microsecond to millisecond range) of such sensors enables one for fast data refreshing. The SAW device itself is small (typical package size 10×5×1 mm 3 ) but the associated radio-frequency (RF) antennas penalize the compactness of the whole sensor, depending on its frequency operation and the nature of its environment.

For instance, operation in dense absorbing (organic) media or with metallic surroundings still is a challenge for achieving large interrogation distance.

We have here investigated the use of SAW resonators and delay lines as buried sensors for long term temperature monitoring: while interrogation speed is hardly an issue, the interrogation distance will define the system efficiency and the range of use. For interrogation distances smaller than 1 meter, applications mainly concern concrete surface properties monitoring [START_REF] Clemena | Handbook on nondestructive testing of concrete[END_REF] , road aging or near surface soil properties monitoring. However, the application range is greatly enhanced if tens-of-meters-range interrogation distances can be reached [START_REF] Gogineni | [END_REF] , since deep soil properties then can be accessed. We report on the long-term monitoring of soil temperature using the two above-mentioned sensor configurations. Although we demonstrate long term stability and robustness of packaged sensors and signal-to-noise ratio compatible with the expected application, the interrogation range (max 80 cm) is insufficient for most geology or geophysics purposes. We then draw our inspiration from the literature concerning Ground Penetrating Radar (GPR) [START_REF] Daniels | Ground Penetrating Radar[END_REF] . The latter technique is widely used for monitoring dielectric interfaces in buried structures, with penetration distances depending on the probe electromagnetic pulse duration and dielectric properties of the soil. We focus on providing complementary informations from sensors with interrogation techniques compatible with GPR, following a strategy commonly known as cooperative target 10 . We particularly focus on the use of delay lines, as the corresponding interrogation method is similar to the one used by GPR which allows for interrogation distances ranging from a few meters to tens of meters and which operates in the lower RF range, depending on soil water content, permittivity and conductivity. Assuming propagation losses in a pure dielectric medium with negligible conductivity (snow or ice), an interrogation distance of about 40 m is predicted which reveals compliant with geology and geophysics purposes (temperature and stress monitoring for instance). These results were experimentally tested, using a delay line designed to operate around 100 MHz (the actual GPR working frequency) exhibiting a very simple time response (3 bits) as coding was not a purpose of this work. The delay line has been buried in snow and interrogated at various depths up to the maximum experimentally feasible distance of 5 m, which tends to validate the predicted interrogation distance.

II. BURIED RESONATORS AS PASSIVE TEMPERATURE SENSORS

A first set of experiments has been performed using resonator based sensors [11][12][START_REF] Kalinin | IEEE Ultrasonics Symposium[END_REF] . Three 434-MHz surface acoustic wave sensors were buried in clay after being connected to dipole antennas. The length of these antennas is adjusted prior to installation in soil assuming a relative permittivity of 10. The purpose of this experiment is to validate the operation of sensors buried in soil and the evolution of the RF link quality over time, as a function of temperature or climatic conditions (for instance moisture level in soil).

Each sensor is made of two resonators connected in parallel, one reference frequency and one measurement frequency within the 1.7-MHz wide European ISM band, with each resonator designed so that its frequency remains within one half of the allocated frequency band (for temperatures ranging from -20 o C to 160 o C). Each sensor is packaged and hermetically sealed in 5×5 mm 2 ceramic packages. The gold-coated contact pads are tin-soldered to the antennas made of 1.6 mm thick FR4 epoxy coated with 30 µm copper. Interrogating these sensors is performed using a custom-designed monostatic pulse mode RADAR system acting as a reflection-mode frequency-sweep network analyzer. The pulse mode operation improves the isolation between the emission and reception steps and hence the interrogation distance, typically a few meters in air [START_REF] Friedt | [END_REF] .

The first observation during installation of the experiment is that an interrogation unit generating 10 dBm, with a detection limit of -70 dBm [START_REF] Friedt | [END_REF] , is unable to detect a usable signal from devices buried only 30 cm deep, consistent with the tabulated electromagnetic propagation losses in clay whose relative permittivity ε r is in the [4..40] range and conductivity σ is in the [2 × 10 -3 ..1] S/m range of 1 to 300 dB/m 15 , as computed through the exponentially decay loss α of a monochromatic electromagnetic plane wave at pulsation ω propagating in a conducting medium

α = ε r 2   1 + σ 2 ε 2 ω 2 -1   where ε = ε 0 ε r .
The read out distance is increased by inserting an electromagnetic waveguide (a simple conducting wire) in the hole in the soil near the buried device (Fig. 1). It must be noted that no electrical connection is provided between this metallic wire and the sensor on one side, or the interrogation unit on the other side, meaning that this setup is resistant to soil motion, oxidation or surface disturbances such a lawn mowing (Fig. 1). Another sensor is then buried at a depth of 80 cm and soldered to a RG174 coaxial cable protruding from the hole in the ground as an open feed connexion. All 10 cm diameter holes were refilled with the same clay than the surrounding area and watered to avoid any air gap.

This setup provides relative temperature informations (Fig. 2) over time as the sensors had not been calibrated prior to the experiment. Long term drift due to aging of the transducers is reduced through the differential measurement approach: although a single resonator frequency might drift over time due to surface contamination after packaging [START_REF] Shreve | Ultrasonics Symposium Proceedings[END_REF][START_REF] Parker | Ultrasonics Symposium Proceedings[END_REF] , the differential approach of measuring a frequency difference between two resonators submitted to the same environment reduces this effect.

The evolution of the temperature provided by the buried sensors is consistent with a sliding average over 9 days of surface temperatures as provided on the web site http: //www.meteociel.fr/ (maximum of the cross-correlation between the experimental data and the averaged values as a function of sliding window length). The result of this experiment running for more than 500 days is exhibited in Fig. 2: the SAW sensors packaged in ceramic housings are resistant to environmental corrosion, and no significant drift or signal loss was observed during the experiment period. The error bars are consistent with a sub-kelvin resolution, typically of the order of ±0.1 K. The efficiency of the wireless link was validated during night-time measurements (no visual identification of the location of the sensor other than by scanning the interrogation unit antenna over a ∼1 m 2 area where the sensor was supposed to be located until a usable RF signal was acquired) or when snow was covering the measurement area.

However, due to the high duty cycle of resonator interrogation (typically 50% emission and 50% reception and signal processing), peak and average RF power are both in the tens to hundreds of milliwatt range, reducing the interrogation range if RF emission regulations are met. On the other hand, ultra-wide band pulse mode RADAR exhibits very low duty cycles, typically 0.1%, associated with high peak powers in the hundreds to thousands of watts. As this interrogation mode is hardly compatible with resonator-based sensor operation, we have considered the possibility of using wide-band SAW devices, i.e. delay lines built on lithium niobate, to meet our goal. Hence, we consider in the next section the use of a commercially available GPR unit as interrogation units for buried acoustic delay lines acting as sensors.

III. INTERROGATING DELAY LINES A. GPR operation

Throughout this presentation, we will focus on the read out of a SAW sensor using a Malå Geoscience (Malå, Sweden) RAMAC GPR equipped in a 100 MHz bistatic configuration.

The simplest implementation of RADAR interrogation units are designed to generate a short -ideally single -pulse including as much energy as possible. This result is achieved in the RF range by slowly loading a capacitor with a high voltage (provided by a switching power supply for embedded designs) and "instantaneously" emptying this energy in an antenna through an avalanche transistor when triggered by a clock pulse. The duration of the energy transfer is defined by the antenna impedance, which is itself influenced by the antenna dimensions and surrounding medium permittivity (Fig. 3). Hence, GPR operation should be considered as fixed wavelength (defined by the antenna dipole dimensions) rather than fixed frequency, since the soil permittivity affects the electromagnetic velocity and hence the pulse central frequency.

In a classical mode of operation, the bistatic GPR unit operates as follows:

1. a radiofrequency pulse is generated by the emitter, for example by triggering the base of an avalanche transistor and letting the current flow from a capacitor loaded with a high voltage (360 V in the case of the RAMAC unit) to the emitting antenna. In this particular case, the peak power in the dipole antenna load (70 Ω impedance at resonance) is thus 2 kW. For a given position of the GPR, a series of time domain return signals is called a GPR trace. Presenting the returned signal power (trace) as a color or grey-scale map is called a scan. Displaying multiple scans side by side for various positions of the GPR unit is called a GPR profile. GPR profile usually maps the evolution over distance of a dielectric interface or obstacle, for example a glacier bedrock (Fig. 3). Our GPR unit performed in agreement with results found in the literature [START_REF] Barrett | [END_REF] , allowing for the identification of an usable signal more than 150 m deep when used on ice to monitor the bedrock interface of a glacier (Fig. 3).

B. GPR for probing acoustic delay lines

Any impedance mismatch between the avalanche transistor output and antenna through a balun will induce ringing and, in classical RADAR applications, unwanted additional oscillations beyond the main pulse. This ringing may be suitable for interrogating delay lines since more than a single pulse is necessary to efficiently load energy into SAW devices, as their pass-band rarely overpasses 10% of their central operating frequency (related to their electro-mechanical coupling). The extreme case is the resonator of quality factor Q which needs Q/π periods (at 434 MHz, Q ≃10000 which yields about 3500 periods) to be efficiently loaded. The quality factor of the antenna is usually much below this value, of the order of unity, and hence a passive resonator (coaxial line) might be added between the balun and the antenna to store energy and induce enough ringing when interrogating resonators. Furthermore, it is wise to detune the antenna and the resonator to avoid too strong a coupling between these elements 19 .

Hence, interrogating delay lines using a RADAR setup includes new challenges. The number of oscillations of the emitted pulse as well as the central frequency are strongly dependent on the permittivity of the surrounding medium. Indeed, the fixed quantity is the emitted signal wavelength which is defined by the size of the dipole antenna of the GPR, while the center frequency is induced by the equivalent permittivity of the air-soil interface. Despite the impact of the environment on the emitted signal, we observe that the emitted signal is so broad in the frequency domain that it will always overlay the relatively narrowband response of the delay line (Fig. 4).

The sensor we have designed includes a transducer made of 21 IDT pairs, three mirrors also made of 21 IDT pairs located at distances from the transducer so that the reflected echoes are detected 1.0, 1.3 and 2.8 µs after the excitation pulse. The acoustic wavelength of λ = 40 µm yields a central frequency around 100 MHz, matching the pulse length generated by the 100 MHz antenna of the RAMAC GPR unit. The 128-degrees Y-rotated black lithium niobate (pyro-free) substrate was selected for its strong piezoelectric coupling as well as large temperature drift, making it ideal for temperature measurement applications. The free surface acoustic velocity of the Rayleigh wave is 3979 m/s: the delay line aluminum grating parameters correspond to a metallisation ratio a/p = 0.5 and a relative height h/λ=2.5%

(1 µm thick aluminum layer) (Fig. 4, top-right) 20 . Although the acoustic sensor itself is less than 1×1 cm 2 in dimensions, the associated 100 MHz antenna is made of a 1 mm-diameter copper-wire dipole of total length 75 cm.

Furthermore, the time stretching strategy used by GPR to achieve such high sampling rates with rather basic electronics is interesting to develop 21 : successive pulses are generated and the response of the environment is monitored after a programmable time delay, which, in the case of a 500 MHz sampling rate, is increased by 2 ns steps at each interrogation iterate. In the case of Malå's RAMAC GPR, the repetition rate is 100 kHz: this time delay of 10 µs between each emitted pulse explains that under favorable conditions, some leftover echo signal from the delay line (up to 10 µs after the excitation pulse has been received by the sensor) is visible before the excitation pulse is emitted. This repetition rate also defines the maximum time delay of the last echo generated by the delay line sensor. Most interesting to our signal processing strategy, this measurement technique allows for fast sampling at baseband of the received signal as opposed to a demodulated magnitude information which might lack the phase information we will be using for determining accurately the time delay between the emitted excitation pulse and received echoes.

Figs. 5 and6 show, on the top-left chart, the time evolution of the reflected signal for a sensor located at 50 cm (Fig. 5) and 1 m (Fig. 6) from the receiving antenna: the sensor is located on the surface of a concrete area, away from the emitting antenna. These sensors were heated up from room temperature to 80 o C by a 3 Ω power resistor supplied with a 1.5 A current: the temperature T was monitored using a Pt100 temperature probe glued next to the acoustic delay line. Since the delay lines were patterned on a (YXl)/128 o lithium niobate cut with an experimental temperature drift coefficient of -70±2 ppm/K, the echo delay ∆τ variation is ∆τ = 70 × 10 -6 τ (T -T 0 ) with T 0 the reference temperature, τ = 1.0 or 1.3 µs depending on the reflection delay τ under consideration. This time delay, ∆τ ∈ [5.5..6.7] ns, is observed as a magnitude signal shift of 3 pixels at most when sampled at 500 MHz. However, as shown by Reindl & al. [22][23][START_REF] Stelzer | Second International Symposium on Acoustic Wave Devices for Future Mobile Communication Systems[END_REF] , the magnitude information provides a rough estimation of the temperature while the use of the phase in that purpose improves the accuracy, although with a modulo 2π uncertainty. One full phase rotation is easily identified in practical conditions: considering τ ≤ 2 µs, a phase rotation of 2π occurs, in our case, when ∆τ =10 ns, which happens when ∆T = 71 K.

We have thus applied the following algorithm to extract the temperature information from the RADAR recordings:

• roughly identify the echo location using a cross-correlation magnitude maximum between the emitted pulse and the received echoes. The first three maxima are considered since we know our delay line is designed with three reflectors, while four echoes are actually seen in the 0-3.5 µs time interval due to additional reflections on the edges of the chip,

• perform the Fourier transform of the returned echo to identify the frequency range of interest

• the accurate time delay deduced from the position of the whole echo burst position is accessible through the phase of the short-term Fourier transform. This value is plotted in the bottom-left graphs of Figs. 5 and6, and compared to the Pt100 temperature probe recording.

The absolute phase of the Fourier transform, i.e. absolute position of the echoes, depends on the distance of the sensor to the receiving antenna, as can be seen at trace 400 of Fig. 6, where the sensor was moved from 1 m from the receiving antenna to 50 cm. Not only does the magnitude of the received signal increase, but more significantly the phase shift of both echoes is affected by this signal change. Hence, the phase difference between the time delays of the two echoes due to the two mirrors on a same delay line yields a reliable estimate of the temperature variations, and absolute temperature when calibrated. However, since the noise of the two time delay estimates are uncorrelated, the noise level of the phase difference is equal to the sum of the noises of each phase estimate: while each phase measurement allows to estimate a temperature with sub-kelvin accuracy when the sensor is located at a fixed position, the temperature recorded from a phase difference is accurate to a standard deviation of 1.5 K when the sensor is located at a distance of 50 cm from the receiving antenna. This figure degrades when the sensor is moved away from the receiving antenna, to get above 3 K when the sensor is located at 1 m from the receiving antenna (Fig. 6, bottom left, red). This short term noise is strongly reduced by stacking multiple estimates in a sliding average, as can be seen in Fig. 6 where the green line is a sliding average over 10 samples of the red phase-to-temperature conversion from single measurements.

These experiments were performed on concrete, with the distance between the emitting antenna and the receiving antenna equal to 1 m, and the sensor located 50 cm or 1 m from the receiving antenna, away from the emitting antenna (which was thus located 1.5 to 2 m from the delay line sensor). Such a configuration is not favorable for efficient coupling since the surface electromagnetic wave is weak with respect to the electromagnetic coupling towards the soil thanks to its strong dielectric permittivity.

Hence, experimenting in a condition favorable to GPR with an environment of low conductivity, provides convenient conditions to assess the practical usage range of these sensors.

We buried sensors in 5 m high snowdrifts close to Ny-Alesund (Spitsbergen, Norway) as a representative environment of temperature monitoring of a glacier in a polar environment.

The signal to noise ratios up to this depth allows for extracting the echoes returned from the SAW delay line, identifying the relative phase values and hence the temperature (Fig. 7).

D. Range estimate and sensor signal identification

The unusually long delay between the incoming pulse and the echoes returned by the delay line -1 to 2 microseconds would be associated with reflectors 85 to 170 m deep in ice -allows in most situations for time domain multiplexing, with sensor-associated signals observed in a time window inconsistent with dielectric interfaces. This method for identifying the source of the signal -dielectric interface or acoustic sensor -is reminiscent of Time Division Multiple Access (TDMA) classically user for sharing a single transmission canal amongst multiple applications.

Considering the usable reflections recorded from ice-rock interfaces more than d interf ace =100 m below the surface (Fig. 3), we wish to estimate the depth [START_REF] Leucci | Scholarly research exchange[END_REF] at which a GPR-like interrogation scheme would be able to detect informations from a buried delay line. Based on the reflection coefficient of the permittivity mismatch at the interface between the two layers and the typical insertion loss of delay lines, we can estimate the range at which a delay line will provide the receiver of the RADAR with enough power for a measurement:

• assuming a plane wave reaching an interface between ice and rock, the Fresnel reflection coefficient R is computed using relative permittivities ε ice = 3.1 [START_REF] Jiang | [END_REF] and ε rock ≃ 5 as

R = √ ε ice - √ ε rock √ ε ice + √ ε rock 2
. We deduce that in this case, the ice-rock interface exhibits an

IL interf ace = 19 dB reflection coefficient
• the ice-rock interface hence presents a reflection coefficient much larger than the typical delay line with a S 11 insertion loss at 35 dB 5 (Fig. 4), meaning that the delay line must be close to the RADAR to provide a meaningful signal

• the Free Space Propagation Loss (FSPL 27 ) calculation is adapted to the RADAR configuration considering that the SAW target acts as a point like source. Hence, the classical Friis formula stating that the electromagnetic power decays as the distance d squared becomes a fourth power law, while the antenna aperture remains proportional to the electromagnetic wavelength λ: F SP L RADAR (d) = 10 × log 10

λ 2 4π × 1 (4πd 2 ) 2 = 10 log 10 λ 2 (4π) 3 d 4 .
• in order to assess the range at which the delay line with IL SAW = 35 dB insertion loss can be interrogated, we must identify the depth d SAW for which the received power is equal to the one computed previously in the case of the reflection on the bedrock:

F SP L RADAR (d interf ace )+IL interf ace = F SP L RADAR (d SAW )+IL SAW yielding
a computation of d SAW independent on λ and numerical constants:

d SAW = d interf ace × 10 (IL interf ace -IL SAW )/40
In our case, since IL SAW -IL interf ace =16 dB, we conclude that the depth at which the acoustic delay line echoes are of the same magnitude than the reflected signal from an ice-rock interface is d SAW ≃ 40 m.

The conclusion of this plane wave analysis is that a SAW delay line buried in ice at a depth of 40 m should provide the same signal level than the dielectric interface at 100 m.

The delay line signature in an echo v.s. antenna position graphics (as shown in Fig. 7 for example) is characterized by multiple hyperbolas translated in time towards greater depths since the acoustic signal is an attenuated replica of the electromagnetic pulse delayed a few microseconds in time. An inter-correlation between the various pulses thus allows an accurate identification of the time delays within the delay line and hence identification of the physical quantity affecting these delays.

Locating the sensor position during a GPR scan is possible through the identification of the hyperbola summit: the reflected signal delay is minimum when the antennas are positioned above the sensor. However, considering a homogeneous medium (in our case ice) with a known electromagnetic velocity c, then the hyperbola equation of the two way time travel 2t as a function of antenna position x on the surface is

4c 2 t 2 -x 2 = d 2
for a sensor located at depth d. Hence, beyond the spatial position of the sensor obtained by scanning the GPR instrument, the depth of the point-like sensor is indicated by the hyperbola curvature equal to 1 d×c . Furthermore, this curvature provides a unique signature response since the observed delay (including the acoustic delay of several hundreds of nanoseconds, which would account for a depth of several tens of meters if it were due to the electromagnetic propagation speed) is inconsistent with a dielectric reflector located at depth d.

Finally, multiple sensors with different polarizations can be located in common view of the GPR unit 18 : we have observed that the strong linear polarization of the pulse emitted by the GPR dipole is able to select the response from a single sensor buried 2 m deep in snow without interference from another similar sensor located about 4 m away and positioned with an orthogonal polarization. This strategy is only possible in the far field range, at a distance of several wavelengths (1.7 m at 100 MHz in ice) from the surface.

IV. CONCLUSION

We have demonstrated that SAW resonators packaged in ceramic packages buried in clay can operate for more than one year with no significant drift or signal quality degradation. Systematic monitoring of these buried devices provides temperature evolutions consistent with surface temperatures. We also have shown that an interrogation unit compliant with the 434 MHz European ISM band allows for interrogating buried sensors at a depth of 80 cm using a coaxial connection, and of 60 cm by promoting the electromagnetic field penetration in soil using a simple conductive wire placed near the sensor and standing out the ground.

We have designed and fabricated a dedicated temperature sensor for use with a 100 MHz ground penetrating RADAR (GPR) unit and demonstrated the ability to record echo signals when the sensor is located at the surface 50 cm and 1 m away from the receiving antenna, as well as buried more than 5 m deep in snow or ice. We have developed the signal processing steps from raw GPR data to extract a temperature informations deduced from the time relative time delay between successive echo pulses.

In order to improve the interrogation depth of sensors, we have analyzed the interrogation strategy of ground penetrating RADARs (GPR), able to detect informations of reflected electromagnetic energy at dielectric interfaces up to 100 m deep at 100 MHz in low loss propagation media such as ice. We extend this result to an estimate of the depth at which a SAW delay line might provide the same amount of reflected energy by compensating the large insertion loss by bringing the sensor closer to the surface: a plane wave calculation of Fresnel reflection coefficient hints a possible depth of 40 m, in agreement with the observed signal to noise ratio achieved when the sensor is located 5 m under the surface. the depth at which the minimum signal will be detectable. Considering a signal to noise ratio above 1.7, the maximum readout distance should be 8 times further than the current position (since 1.7 4 ≃ 8), consistent with the 40 m maximum depth estimated from the classical GPR link budget presented in the main text.

2 .

 2 the direct electromagnetic wave propagating on the surface, as well as all the echoes reflected from the dielectric interfaces in the ground, are recorded by the receiving unit at a sampling rate at least 10 times the nominal value of the emitted pulse (in this case a nominal working frequency of 100 MHz), with a sampling triggered by the same signal controlling the base of the avalanche transistor 3. equivalent time sampling reduces the receiving unit cost and bandwidth: the emitted pulse is repeated at a rate slower than the inverse time needed for the pulse to reach the maximum probing depth (10 µs repetition rate in the case of the RAMAC unit, yielding a maximum probing depth of 850 m considering an electromagnetic velocity in ice of 170 m/µs) and the returned signal is recorded after a time interval referenced on the trigger signal and increased by time steps inverse of the wanted sampling frequency. Hence, by delaying the recording time by an additional 250 ps with respect to the trigger signal every new emitted pulse, an equivalent sampling rate of 4 GHz is achieved even with much slower analog to digital converters and low communication bandwidth between the receiving unit and the recording computer.

FIG. 1 :FIG. 2 :FIG. 3 :FIG. 4 :

 1234 Figures

FIG. 5 :

 5 FIG.5: Top left: the raw color-coded time evolution of the recorded RADAR echo magnitude between 1.0 µs and 1.3 µs after the excitation pulse was emitted. The sampling is performed at 500 MHz, or 5 times the frequency of the signal of interest. Top-right: identification of the frequency component (index) representative of the delay line, here visible as a maximum of the magnitude of the Fourier transform of the points from 0.9 to 1.1 µs (first echo) and 1.2 to 1.4 µs (second echo). We observe that this frequency component of interest does not change with temperature (i.e. is independent on the trace number). Bottom right: time-evolution of the unwrapped phase of the Fourier transform at frequency abscissa 25 as identified from the top-right graph. Bottomleft: time evolution of the phase difference between the first and second echoes, after scaling and translation to match the reference temperature curve recorded with a Pt100 probe located next to the delay line. During this whole experiment, the receiving antenna is located 1 m from the emitting antenna, and the sensor is 50 cm from the receiving antenna away from the emitting antenna.

FIG. 6 :

 6 FIG.6:The graph sequence and analysis is the same than the one described in the caption of Fig.5. Here, however, the sensor is first located 1 m from the receiving antenna, away from the emitting antenna, and brought closer to 50 cm of the receiving antenna at trace number 400. This distance change is observed as an increase of the magnitude of the signal of interest (top right graph, magnitude of the Fourier transform of the echo), a phase shift in the bottom graph affecting both echoes in the same way, and a decrease of the temperature estimate standard deviation (bottom-left graph).

FIG. 7 :

 7 FIG.7: Signal acquired while scanning a 100 MHz GPR unit over a sensor buried 2.20 m deep in snow. The emitted pulse exhibits ringing due to impedance mismatch, a condition degrading depth resolution but favorable to efficiently load the acoustic delay line. The recorded signal clearly displays four echoes, the first three being used to extract the physical quantity under investigation. The absolute phase with respect to the emitted pulse is dependent on antenna position and constantly rises as the RADAR is brought close to the sensor, but the phase difference is independent on antenna position and is representative of the physical quantity under investigation. As expected, each echo is made of 21 oscillations, which is equal to the number of electrode pairs in the transducer. The "inverted" hyperbola shape of the echoes is an aliasing artifact when displaying the data.

FIG. 8 : 7 FIG. 9 :

 879 FIG.8: Experimental setup for recording signals from a sensor while scanning a 100 MHz GPR unit over a sensor buried 5 m deep in snow. The emitted pulse exhibits ringing due to impedance mismatch, a condition degrading depth resolution but favorable to efficiently load the acoustic delay line. The recorded signal clearly displays four echoes, the first three being used to extract the physical quantity under investigation.

Acknowledgment J.-M Friedt would like to thank S. Zhgoon (Moscow Power Engineering Institute, Russia) for fruitful discussions during the 2009 IFCS/EFTF conference. Ground Penetrating Radar measurements on glaciers were performed as part of French National Research Agency (ANR) Hydro-Sensor-FLOWS program under the supervision of D. Laffly, C. Marlin, M.

Supplementary material

Fig. 8 depicts the experimental setup for measuring the response of SAW delay lines buried 5 m deep in a snowdrift. The sensor and the associated 70 cm long dipole antenna are located in a 3 cm diameter tube filled with snow. The tube is inserted about 1.5 m deep in snow while the GPR scans this area and records both dielectric reflections and echoes from the SAW sensor as a function of antenna position. Although each absolute echo phase with respect to the emitted pulse is dependent on the antenna position, the difference of the phases of the echoes is independent on the antenna position and only depends on the acoustic velocity, or in this case the temperature through the temperature coefficient of frequency of the piezoelectric substrate. The data displayed in Figs. 7 and8 were processed using the Seismic Unix package (http://www.cwp.mines.edu/cwpcodes) with the application of a normalization step and bandpass filtering in the 100±50 MHz band.

Analyzing the signal to noise ratio of the echoes returned -starting 1 µs after the emitted pulse and for a duration of 3.3 µs -from the sensor buried 5 m deep in snow, one can estimate