Didier Lime

Olivier H Roux

Charlotte Seidner

Parameter Synthesis for Bounded Cost Reachability in Time Petri Nets

We investigate the problem of parameter synthesis for time Petri nets with a cost variable that evolves both continuously with time, and discretely when firing transitions. More precisely, parameters are rational symbolic constants used for time constraints on the firing of transitions and we want to synthesise all their values such that the cost variable stays within a given budget. We first prove that the mere existence of such values for the parameters is undecidable. We nonetheless provide a symbolic semi-algorithm that is proved both sound and complete when it terminates. We also show how to modify it for the case when parameters values are integers. Finally, we prove that this modified version terminates if parameters are bounded. While this is to be expected since there are now only a finite number of possible parameter values, this is interesting because the computation is symbolic and thus avoids an explicit enumeration of all those values. Furthermore, the result is a symbolic constraint representing a finite union of convex polyhedra that is easily amenable to further analysis through linear programming. We finally report on the implementation of the approach in Romeo, a software tool for the analysis of hybrid extensions of time Petri nets.

Introduction

So-called priced or cost timed models are suitable for representing real-time systems whose behaviour is constrained by some resource consuming (be it energy or CPU time, for instance) and for which we need to assess the total cost accumulated during their execution. Such models can even describe whether the evolution of the cost during the run is caused by staying in a given state (continuous cost) or by performing a given action (discrete cost). Thus, the task of finding if the model can reach some "good" states while keeping the overall cost under a given bound (or, further, finding the minimum cost) can prove of interest in many real-life applications, such as optimal scheduling or production line planning.

Timed models, however, require a thorough knowledge of the system for their analysis and are thus difficult to build in the early design stages, when the system is not fully identified. Even when all timing constraints are known, the whole design process must often be caried out afresh, whenever the environment changes. To obtain such valuable characteristics as flexibility and robustness, the designer may want to relax constraints on some specifications by allowing them a wider range of values. To this end, parametric reasoning is particularly relevant for timed models, since it allows designers to use parameters instead of definite timing values.

We therefore propose to tackle the definition and analysis of models that support both (linear) cost functions and timed parameters.

Related work Parametric timed automata (PTA) [START_REF] Alur | Parametric real-time reasoning[END_REF] extend timed automata [START_REF] Alur | A Theory of Timed Automata[END_REF] to overcome the limits of checking the correctness of the systems with respect to definite timing constraints. The reachability-emptiness problem, which tests whether there exists a parameter valuation such that the automaton has an accepting run, is fundamental to any verification process but is undecidable [START_REF] Alur | Parametric real-time reasoning[END_REF]. L/U automata [START_REF] Hune | Linear parametric model checking of timed automata[END_REF] use each parameter either as a lower bound or as an upper bound on clocks. The reachability-emptiness problem is decidable for this model, but the state-space exploration, which would allow for explicit synthesis of all the suitable parameter valuations, still might not terminate [START_REF] Jovanović | Integer parameter synthesis for real-time systems[END_REF]. To obtain decidability results, the approach described in [START_REF] Jovanović | Integer parameter synthesis for real-time systems[END_REF] does not rely on syntactical restrictions on guards and invariants, but rather on restricting the parameter values to bounded integers. From a practical point of view, this subclass of PTA is not that restrictive, since the time constraints of timed automata are usually expressed as natural (or perhaps rational) numbers.

In [START_REF] Alur | Optimal paths in weighted timed automata[END_REF], the authors have proved the decidability of the optimal-cost problem for Priced Timed Automata with non-negative costs. In [START_REF] Behrmann | Minimum-cost reachability for priced timed automata[END_REF][START_REF] Behrmann | Optimal scheduling using priced timed automata[END_REF][START_REF] Larsen | As cheap as possible: Efficient cost-optimal reachability for priced timed automata[END_REF], the computation of the optimal-cost to reach a goal location is based on a forward exploration of zones extended with linear cost functions. In [START_REF] Bouyer | Symbolic optimal reachability in weighted timed automata[END_REF], the authors have improved this approach, so as to ensure termination of the forward exploration algorithm, even when clocks are not bounded and costs are negative, provided that the automaton has no negative cost cycles. In [START_REF] Abdulla | Priced timed Petri nets[END_REF], the considered model is a timed arc Petri net, under weak firing semantics, extended with rate costs associated with places and firing costs associated with transitions. The computation of the optimal-cost for reaching a goal marking is based on similar techniques to [START_REF] Alur | Optimal paths in weighted timed automata[END_REF]. In [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF], the authors have investigated the optimal-cost reachability problem for time Petri nets where each transition has a firing cost and each marking has a rate cost (represented as a linear rate cost function over markings). To compute the optimal-cost to reach a goal marking, the authors have revisited the state class graph method to include costs.

Our contribution We propose in Section 2 an extension of time Petri nets with costs (both discrete and continuous with time) and timing parameters, i.e., rational symbolic constants used in the constraints on the firing times of transitions.

Within this formalism, we define two problems dealing with parametric reachability within a bounded cost. We prove in Section 3 that the existence of a parameter valuation to reach a given marking under a given bounded cost is undecidable. This proof adapts a 2-counter machine encoding first proposed in [START_REF] Jovanović | Parametric Verification of Timed Systems[END_REF] for PTA. To our knowledge it is the first time a direct Petri net encoding is provided and the adaptation is not trivial. We give in Section 4 a symbolic semi-algorithm that computes all such parameter valuations when it terminates, and we prove its correctness. We propose in Section 5 a variant of this semialgorithm that computes integer parameter valuations and prove in Section 6 its termination provided those parameter valuations are bounded and the cost of each run is uniformly lower-bounded for integer parameter valuations. This technique is symbolic and avoids the explicit enumeration of all possible parameter valuations. The basic underlying idea of using the integer hull operator was first investigated in [START_REF] Jovanović | Integer parameter synthesis for real-time systems[END_REF] for PTA, but this is the first time that it is adapted and proved to work with state classes for time Petri nets, and the fact that it naturally also preserves costs for integer parameter valuations is new and very interesting. We finally describe in Section 7 the implementation of the approach in the tool Romeo by analysing a small scheduling case-study.

Parametric Cost Time Petri Nets

Preliminaries

We denote the set of natural numbers (including 0) by N, the set of integers by Z, the set of rational numbers by Q and the set of real numbers by R. We note Q ≥0 (resp. R ≥0 ) the set of non-negative rational (resp. real) numbers. For n ∈ N, we let 0, n denote the set {i ∈ N | i ≤ n}. For a finite set X, we denote its size by |X|.

Given a set X, we denote by I(X), the set of non empty real intervals that have their finite end-points in X. For I ∈ I(X), I denotes its left end-point if I is left-bounded and -∞ otherwise. Similarly, I denotes the right end-point if I is right-bounded and ∞ otherwise. We say that an interval I is non-negative if I ⊆ R ≥0 . Moreover, for any d ∈ R ≥0 and any non-negative interval I, we let I d be the interval defined by {θ -d | θ ∈ I ∧ θ -d ≥ 0}. Note that this is again a non-negative interval.

Given sets V and X, a V -valuation (or simply valuation when V is clear from the context) of X is a mapping from X to V . We denote by V X the set of V -valuations of X. When X is finite, given an arbitrary fixed order on X, we often equivalently consider V -valuations as vectors of V |X| . Given a V -valuation v of X and Y ⊆ X, we denote by v |Y the projection of v on Y , i.e., the valuation on Y such that ∀x ∈ Y, v |Y (x) = v(x).

Time Petri Nets with Costs and Parameters

Definition 1 (Parametric Cost Time Petri Net (pcTPN)). A Parametric Cost Time Petri Net (pcTPN) is a tuple N = (P, T, P, • ., . • , m 0 , I s , cost t , cost m ) where -P is a finite non-empty set of places, -T is a finite set of transitions such that T ∩ P = ∅, -P is a finite set of parameters, -• . : T → N P is the backward incidence mapping, -. • : T → N P is the forward incidence mapping, m 0 ∈ N P is the initial marking, -I s : T → I(N ∪ P) is the (parametric) static firing interval function, -cost t : T → Z is the discrete cost function, and cost m : N P → Z is the cost rate function.

Given a parameterized object x (be it a pcTPN, a function, an expression, etc.), and a Q-valuation v of parameters, we denote by v(x) the corresponding non-parameterized object, in which each parameter a has been replaced by the value v(a).

A marking is an N-valuation of P . For a marking m ∈ N P , m(p) represents a number of tokens in place p. A transition t ∈ T is said to be enabled by a given marking m ∈ N P if for all places p, m(p) ≥ • t(p). We also write m ≥ • t. We denote by en(m) the set of transitions that are enabled by the marking m:

en(m) = {t ∈ T | m ≥ • t}.
Firing an enabled transition t from marking m leads to a new marking m = m -• t + t • . A transition t ∈ T is said to be newly enabled by the firing of a transition t from a given marking m ∈ N P if it is enabled by the new marking but not by m-• t (or it is itself fired). We denote by newen(m, t) the set of transitions that are newly enabled by the firing of t from the marking m: newen(m, t) = t ∈ en(m

-• t + t • ) | t ∈ en(m -• t) or t = t A state of the net N is a tuple (m, I, c, v) in N P × I(R ≥0 ) T × R × Q P ≥0
, where: m is a marking of N , I is called the interval function and associates a temporal interval to each transition enabled by m. Value c is the cost associated with that state and valuation v assigns a rational value to each parameter for the state.

Definition 2 (Semantics of a pcTPN). The semantics of a pcTPN is a timed transition system (Q, Q 0 , →) where:

-Q ⊆ N P × I(R ≥0 ) T × R × Q P ≥0 -Q 0 = {(m 0 , I 0 , 0, v)|v ∈ Q P ≥0 , ∀t ∈ T, v(I s (t)) = ∅} where ∀t ∈ en(m 0 ), I 0 (t) = I s (t) -→ consists of two types of transitions: • discrete transitions: (m, I, c, v) t∈T --→ (m , I , c , v) iff * m ≥ • t, m = m -• t + t • and v(I(t)) = 0, * ∀t ∈ en(m ) • I (t ) = I s (t ) if t ∈ newen(m, t), • I (t ) = I(t ) otherwise * c = c + cost t (t) • time transitions: (m, I, c, v) d∈R ≥0 ----→ (m, I d, c , v), iff ∀t ∈ en(m), (I d)(t) ≥ 0 and c = c + cost m (m) * d.
A run of a pcTPN N is a (finite or infinite) sequence q 0 a 0 q 1 a 1 q 2 a 2 • • • such that q 0 ∈ Q 0 , for all i > 0, q i ∈ Q, a i ∈ T ∪ R ≥0 and q i ai -→ q i+1 . The set of runs of N is denoted by Runs(N ). We note (m, I, c, v) t@d --→ (m , I , c , v) for the sequence of elapsing d ≥ 0 followed by the firing of the transition t. We denote by sequence(ρ) the projection of the run ρ over T : for a run ρ = q 0 t0@d0 ----→ q 1 t1@d1 ----→

q 2 t2@d2 ----→ q 3 t3@d3 ----→ • • • , we have sequence(ρ) = t 0 t 1 t 2 t 3 • • • . We write q t - → q if there exists d ≥ 0 such that q t@d --→ q .
For a finite run ρ we denote by last(ρ) the last state of ρ and by lastm(ρ) its marking. A state (m, I, c, v) is said to be reachable if there exists a finite run ρ of the net, with last(ρ) = (m, I, c, v). A marking m is reachable for parameter valuation v, if there exists some I and c such that (m, I, c, v) is reachable.

For k ∈ N and parameter valuation v, the (Cost) Time Petri net v(N ) is said to be k-bounded if for all reachable markings m, and all places p, m(p) ≤ k. We say that v(N ) is bounded if there exists k such that it is k-bounded.

The cost cost(ρ) of a finite run ρ, with last state (m, I, c, v) is c. Since we are interested in minimising the cost, the cost of a sequence of transitions σ is defined as cost(σ) = inf ρ∈Runs(N ),sequence(ρ)=σ cost(ρ). For the sake of the clarity of the presentation, we consider only closed intervals (or right-open to infinity) so this infimum is actually a minimum.

Parametric Cost Problems

Given a set of target markings Goal, the problems we are interested in are:

1. the existential problem: Given a finite maximum cost value c max , is there a parameter valuation v such that some marking in Goal is reachable with a cost less than c max in v(N )? 2. the synthesis problem: Given a finite maximum cost value c max , compute all the parameter valuations v such that some marking in Goal is reachable with a cost less than c max in v(N ).

We prove in Section 3 that the existential problem is undecidable.

Undecidability Results

The existential parametric time bounded reachability problem for bounded parametric time Petri nets asks whether a given target marking is reachable for some valuation of the parameter(s) within c max time units. This is a special case of the existential cost bounded reachability problem defined in Section 2, with no discrete cost and a uniform cost rate of 1. Proposition 1 therefore implies the undecidability of that more general problem. Proof. Given a bounded parametric time Petri net N , we want to decide whether there exists some parameter valuation v such that some given marking can be reached within c max time units in v(N ). The idea of this proof was first sketched in [START_REF] Jovanović | Parametric Verification of Timed Systems[END_REF] for parametric timed automata. We encode the halting problem for twocounter machines, which is undecidable [START_REF] Minsky | Computation: Finite and Infinite Machines[END_REF], into the existential problem for parametric timed Petri nets. Recall that a 2-counter machine M has two nonnegative counters (here C x and C y ), a finite number of states and a finite number of transitions, which can be of the form: 1) when in state s i , increment a counter and go to s j ; 2) when in state s i , decrement a counter and go to s j ; 3) when in state s i , if a counter is null then go to s j , otherwise block. The machine starts in state s 0 and halts when it reaches a particular state s halt . Given such a machine M, we now provide an encoding as a parametric time Petri net N M : each state s i of the machine is encoded as place, which we also call s i . The encoding of the 2-counter machine M is as follows: it uses two rationalvalued parameter a and b, and three gadgets shown in Figure 1.a modelling three clocks x, y, z. Recall that, for a state (m, I, c, v), the enabling time of an enabled transition t is v(I s (t) -I(t)). For the gadget modelling the clock x, the value of the clock x is equal to: i) the enabling time of the transition t x=b when P x≤b is marked; ii) b when P x≤b is marked; iii) the sum of b and the enabling time of the transition t x=a+b when P x≥b is marked (note that this value is lower than a + b); iv) a + b when P x=a+b is marked; v) an unknown (an irrelevant) value in all other cases.

[0, 0] tz=b [b, b] εz b [0, 0] tz=a+b [a, a] εz a+b [0, 0] tx=a+b [a, a] tx=b [b, b] ty=a+b [a, a] ty=b [b, b] εx a+b [0, 0] εx b [0, 0] εy a+b [0, 0] εy b [0, 0] start [0, 0] R(y) [0, 0] R(x) [0, 0] done [0, 0] P0 Ploop a P1 Pa Px=b Py=b Pz=0 s0 loopa [a, a] Ploop b Pb loopb [b, b] ta [a, a] tb [b, b] t1 [1, 1] εa [0, 0] εb [0, 0] ε1 [0, 0] start [0, 0]
The gadget encoding the increment instruction of C y is given in Figure 1.a. The clocks x and y store the value of each counter C x and C y as follows x = b -a.C x and y = b -a.C y when z = 0. The zero-test gadget is given in Figure 2. We use the initial gadget in Figure 1.b to initialise a and b such that 0 < a ≤ 1 and 0 < b ≤ 1. The system is studied over 1 time unit.

Increment: We start from some encoding configuration: x = b -a.C x , y = b -a.C y and z = 0 in a marking such that the places P z=0 and s i are marked. After the firing of the transition start, there is an interleaving of the transitions R(x) and R(y) that go through the gadget. Finally, we can fire the transition done when z = b (i.e. b -a.C x later) and we have z = 0, x = b -a.C x and y = b -a(C y + 1) as expected. Moreover, v(N M ) will block for all the parameter valuations v which not correctly encode the machine.

Decrement: By replacing the arc from P z=b to done by an arc from P z=a+b to done, the only difference in the previous reasoning is that the elapsing time to fire done is increased of a. Then we obtain z = 0,

x = b+a-a.C x = b-a.(C x -1) and y = b -a.C y corresponding to the decrement of C x .
We can obtain symmetrically (by swapping x and y) the increment of C x and the decrement of C y Both the increment gadget and the zero-test gadget require b time units, and the decrement gadget requires (a + b) time units. Since the system executes over 1 time unit, for any value of a > 0 and b > 0, the number of operations that the machine can perform is finite. We consider two cases: 4 A Symbolic Semi-algorithm for Parameter Synthesis

State Classes

We now introduce the notion of state classes for pcTPNs. It was originally introduced for time Petri nets in [START_REF] Berthomieu | An enumerative approach for analyzing time petri nets[END_REF][START_REF] Berthomieu | Modeling and verification of time dependent systems using time Petri nets[END_REF], and extended for timing parameters in [START_REF] Traonouez | Parametric model-checking of stopwatch Petri nets[END_REF], and for costs in [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF]. We show that those two extensions seamlessly blend together.

For an arbitrary sequence of transitions σ = t 1 . . . t n ∈ T * , let C σ be the set of all states that can be reached by the sequence σ from any initial state q 0 :

C σ = {q ∈ Q|q 0 t1 -→ q 1 • • • tn -→ q}.
[0, 0] tz=b [b, b] εz b [0, 0] tz=a+b [a, a] εz a+b [0, 0] tx=a+b [a, a] tx=b [b, b] ty=a+b [a, a] ty=b [b, b] εx a+b [0, 0] εx b [0, 0] εy a+b [0, 0] εy b [0, 0] start [0, 0] R(y) [0, 0] R(x) [0, 0] 0-test [0, 0]
Fig. 2. Encoding 0-test over bounded-time: when in state si, if Cx = 0 then go to sj and, if we note en(m) = {t 1 , . . . , t n }, then D is the set of points (θ 1 , . . . , θ n , c, v) such that (m, I, c, v) ∈ C σ and for all t i ∈ en(m), θ i ∈ I(t i ). For short, we will often write ( θ, c, v) for such a point, with θ = (θ 1 , . . . , θ n ) and a small abuse of notation. We denote by Θ the set of θ i variables, of which we have one per transition of the net: for the sake of simplicity, we will usually use the same index to denote for instance that θ i corresponds to transition t i .

C σ is called a state class and D is its firing domain. Lemma 1 equivalently characterises state classes, as a straigthforward reformulation of the definition: From Lemma 1, we can then deduce a characterisation of the "next" class, obtained by firing a firable transition from some other class. This is expressed by Lemma 2.

Lemma 2. Let C σ = (m, D) and C σ.t f = (m , D ), we have: ( θ , c , v) ∈ D iff ∃( θ, c, v) ∈ D s.t.        ∀t i ∈ en(m), θ i -θ f ≥ 0 ∀t i ∈ en(m -• t f ), θ i = θ i -θ f ∀t i ∈ newen(m, t f ), θ i ∈ v(I s )(t i ) c = c + cost m (m) * θ f + cost t (t f ) Proof. Consider ( θ , c , v) ∈ D .
Then by Lemma 1, there exists a run ρ in v(N ), and I : en(m) → I(Q ≥0 ), such that sequence(ρ ) = σ.t f , (m , I , c ) = last(ρ ), and θ ∈ I . Consider the prefix ρ of ρ such that sequence(ρ) = σ. The last state of ρ can be written (m, I, c, v) for some I and c. We know that t f is fired from (m, I, c, v) so there exists some delay d such that I(t f ) ≤ d and for all other transitions t i enabled by m, I(t i ) ≥ d. Furthermore, c = c -cost m (m) * dcost t (t f ). It follows that there exists a point θ ∈ I with the desired properties.

The other direction is similar.

Note that according to Lemma 2, D is not empty if and only if there exists ( θ, c, v) in D such that for all t i ∈ en(m), θ i ≥ θ f . In that case we say that t f is firable from (m, D) and note t f ∈ firable((m, D)).

From Lemma 2, it follows that C σ.t f can be computed from C σ using Algorithm 1. Note that it is formally the same algorithm as in [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF].

Given a class C and a transition t firable from C, we note Next(C, t) the result of applying Algorithm 1 to C and t. ; since all the operations on domains in Algorithm 1 are polyhedral, all the domains of state classes are also convex polyhedra. Note that only enabled transitions are constrained in the domain of a state class.

Algorithm 1 Successor (m , D ) of (m, D) by firing t f 1: m ← m -• t f + t • f 2: D ← D ∧ i =f,t i ∈en(m) θ f ≤ θi 3: for all ti ∈ en(m -• t f ), i = f , add variable θ i to D , constrained by θi = θ i + θ f 4: add variable c to D , constrained by c = c + θ f * costm(m) + costt(t f
Naturally, we define the cost of state class C σ as cost(C σ ) = cost(σ).

The Synthesis Semi-algorithm

In Algorithm 2, we explore the symbolic state-space in a classic manner. Whenever a goal marking is encountered we collect the parameter valuations that allowed that marking to be reached with a cost less or equal to c max .

The Passed list records the visited symbolic states. Instead of checking new symbolic states for membership, we test a weaker relation denoted by : does there exist a visited state allowing more behaviors with a cheaper cost?

For any state class C = (m, D) and any point ( θ, v) ∈ D |Θ∪P , the optimal cost of ( θ, v) in D is defined by cost D ( θ, v) = inf ( θ,c,v)∈D c.

( θ, v) ∈ D |Θ∪P , cost D ( θ, v) ≤ cost D ( θ, v).
Algorithm 2 Symbolic semi-algorithm computing all parameter valuations such that some markings are reachable with a bounded cost. The following result is a fairly direct consequence of Definition 3:

Lemma 3. Let C σ1 and C σ2 be two state classes such that C σ1 C σ2 . If a transition sequence σ is firable from C σ1 , it is also firable from C σ2 and cost(C σ1.σ ) ≥ cost(C σ2.σ ). Proof. Let C σ1 = (m 1 , D 1 ) and C σ2 = (m 2 , D 2 ). From Definition 3, for any point ( θ, c 1 , v) ∈ D 1 , there exists a point ( θ, c 2 , v) ∈ D 2 ) such that c 2 ≤ c 1 .
This implies that: (i) cost(C σ1 ) ≥ cost(C σ2 ); (ii) if transition t is firable from C σ1 , then it is firable from C σ2 and Next(C σ1 , t) Next(C σ2 , t). And the result follows by a straightforward induction.

While can be checked using standard linear algebra techniques, we can also reduce it to standard inclusion on polyhedra by removing the upper bounds on cost (an operation called cost relaxation) [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF].

Lemma 4. The following invariant holds after each iteration of the while loop in Algorithm 2: for all C σ = (m, D) ∈ Passed,

1. for all prefixes σ of σ, C σ ∈ Passed; 2. if m ∈ Goal then D ∩ (c ≤ c max ) |P ⊆ PolyRes; 3. if t is firable from C σ -either C σ.t ∈ Waiting, -or there exists C ∈ Passed such that C σ.t C .
Proof. We prove this lemma by induction. Before the while loop starts, Passed is empty so the invariant is true. Let us now assume that the invariant holds for all iterations up to the n-th one, with n ≥ 0, and that Waiting = ∅. Let C σ ∈ Waiting be the selected class at line 5; to check whether the invariant still holds at the end of the (n + 1)-th iteration, we only have to test the case where C σ is added to Passed (which means that the condition at line 9 is true). We can then check each part of the invariant:

1. C σ was picked from Waiting (line 5); except for the initial class (for which σ is empty, and therefore has no prefix), it means that, in a previous iteration, there was a sequence σ and a transition t ∈ firable(C σ ) such that σ = σ .t Both the basis case and the induction step are true: the result follows by induction.

(
Proposition 2. After any iteration of the while loop in Algorithm 2:

1. if v ∈ PolyRes, then there exists a run ρ in v(N ) such that cost(ρ) ≤ c max and lastm(ρ) ∈ Goal. 2. if Waiting = ∅ then, for all parameter valuations v such that there exists a run ρ in v(N ) such that cost(ρ) ≤ c max and lastm(ρ) ∈ Goal, we have v ∈ PolyRes. From Lemma 3, sequence σ 3 is also firable from C σ and C σ .σ3 = (m, D ), with cost(C σ .σ3 ) ≤ cost(C σ ) ≤ c max . By Lemma 1, there exists thus a run ρ in v(N ), with sequence(ρ ) = σ .σ 3 , lastm(ρ ) ∈ Goal and cost(ρ ) ≤ c max . Also, from Lemma 4 (item 1), we know that for all prefixes of σ , the corresponding state class is in Passed, so the biggest suffix of σ .σ 3 as defined above in the induction hypothesis has length less or equal to n -1, and the induction hypothesis applies to ρ , which allows to conclude.

In particular, if the algorithm terminates, then the waiting list is empty and PolyRes is exactly the solution to the synthesis problem.

Restricting to Integer Parameters

Obviously, in general, (semi-)Algorithm 2 will not terminate, since the emptiness problem for the set it computes is undecidable.

To ensure termination, we can however follow the methodology of [START_REF] Jovanović | Integer parameter synthesis for real-time systems[END_REF]: we require that parameters are bounded integers and, instead of just enumerating the possible parameter values, we propose a modification of the symbolic state computation to compute these integer parameters symbolically. For this we rely on the notion of integer hull.

We call integer valuation a Z-valuation. Note that a Z-valuation is also an R-valuation, and given a set D of R-valuations, we denote by Ints(D) the set of integer valuations in D.

The convex hull of a set D of valuations, denoted by Conv(D), is the intersection of all the convex sets of valuations that contain D.

The integer hull of a set D of valuations, denoted by IH(D), is defined as the convex hull of the integer valuations in D: IH(D) = Conv(Ints(D)).

For a state class C = (m, D), we write IH(C) for (m, IH(D)).

Before we see how our result can be adapted for the restriction to integer parameter valuations, and from there how we can enforce termination of the symbolic computations when parameters are assumed to be bounded, we need some results on the structure of the polyhedra representing firing domains of cost TPNs.

By the Minkowski-Weyl Theorem (see e.g. [START_REF] Schrijver | Theory of linear and integer programming[END_REF]), every convex polyhedron can be either described as a set of linear inequalities, as seen above, or by a set of generators. More precisely, for the latter: if d is the dimension of polyhedron P , there exists v 1 , . . . , v p , r 1 , . . . , r s ∈ R d , such that for all points x ∈ P , there exists λ 1 , . . . , λ p ∈ R, µ 1 , . . . , µ s ∈ R ≥0 such that i λ i = 1 and x = i λ i v i + i µ i r i . The v i 's are called the vertices of P and the r i 's are the extremal rays of P . The latter correspond to the directions in which the polyhedron is infinite. In our case, they correspond to transitions with a (right-)infinite static interval, and possibly the cost.

A classic property of vertices, which can also be used as a definition, is as follows: v is a vertex of P iff for all non-null vectors x ∈ R d , either v + x ∈ P or v -x ∈ P (or both), + and -being understood component-wise.

Proposition 3. Let N be a (non-parametric) cost TPN and let C = (m, D) be one of its state classes, then D has integer vertices.

Proof. We have proved in [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF] that the domain D of a state class of a cost TPNs, with removed upper bounds on cost (so-called relaxed classes), can be partitioned into a union of simpler polyhedra n i=1 D i that have the following key properties: (1) by projecting the cost out we obtain a convex polyhedron D i|Θ with integer vertices (actually a zone, as in [START_REF] Larsen | Model-Checking for Real-Time Systems[END_REF][START_REF] Berthomieu | Modeling and verification of time dependent systems using time Petri nets[END_REF]), and (2) these simpler polyhedra all have exactly one constraint on the cost variable, i.e., of the form c ≥ ( θ), with integer coefficients. Note that the same result can be obtained, with the same technique, if we consider non-relaxed state classes, except that, we also have an upper bound on cost that is always greater or equal to the lower bound. We prove in Lemma 5 that each of these simpler polyhedra also has integer vertices. Since D and each of the D i 's are convex and since D = i D i , D is equal to the convex hull of the vertices of the D i 's and therefore D also has integer vertices. Lemma 5. Let D be a convex polyhedron on variables θ 1 , . . . , θ n , c such that the projection of P on the θ variables has integer vertices, and there are two constraints on c of the form c ≥ (θ 1 , . . . , θ n ) and c ≤ (θ 1 , . . . , θ n ), with and linear terms with integer coefficients, such that (θ 1 , . . . , θ n ) ≤ (θ 1 , . . . , θ n ), for all values of the θ i 's.

Then, the vertices of D are the points (θ 1 , . . . , θ n , (θ 1 , . . . , θ n )) and (θ 1 , . . . , θ n , (θ 1 , . . . , θ n )) such that (θ 1 , . . . , θ n ) is a vertex of D |Θ , and they are integer points.

Proof. Recall here that we consider all constraints in D to be non-strict so all polyhedra are topologically closed. The reasoning extends with no difficulty to non-necessarily-closed polyhedra by considering so-called closure points in addition to vertices [START_REF] Bagnara | Not necessarily closed polyhedra and the double description method[END_REF].

Consider a non-vertex point θ in D |Θ and let ( θ, c) be a point of D. Then using the form of the unique cost constraint, we have c ≥ ( θ). Now since θ is not a vertex, there exists a vector x such that both θ + x and θ -x belong to D |Θ . Then, for sure, ( θ + x, ( θ + x)) ∈ D and ( θ -x, ( θ -x)) ∈ D. And since is linear, ( θ + x, ( θ) + ( x))) ∈ D, i.e., ( θ, ( θ)) + ( x, ( x)) ∈ D. And similarly, ( θ, ( θ)) -( x, ( x)) ∈ D. Using again the form of the unique cost constraint, and the fact that c ≥ ( θ), we finally have ( θ, c)

+ ( x, ( x)) ∈ D and ( θ, c) -( x, ( x)) ∈ D, that is, ( θ, c) is not a vertex of D.
By contraposition, any vertex of D extends a vertex of D |Θ , and using a last time the form of the cost constraint, any vertex of D, is of the form ( θ, ( θ)), with θ a vertex of D |Θ : suppose ( θ, c) is a vertex of D, with c > ( θ), then for

x defined with c -( θ) on the cost variable, and 0 on all other dimensions, we clearly have both ( θ, c) + x and ( θ, c) -x in D, which is a contradiction.

We conclude by remarking that, since D |Θ has integer vertices, all the coordinates of θ are integers, and since has integer coefficients then ( θ) is an integer.

We can deal with the upper bound defined by in exactly the same way.

From Proposition 3, we can prove the following lemma that will be very useful in the subsequent proofs. When we restrict ourselves to integer parameter but continue to work symbolically, we need to adjust the definitions of the firability of a transition from a class and of the cost of a class.

First, a transition t f is firable for integer parameter valuations from a class (m, D), call this N P -firable, if there exists an integer parameter valuation v and a point ( θ, c, v) in D such that for all transitions t i ∈ en(m), θ i ≥ θ f . Lemma 7. Let C = (m, D) be a state class. Transition t f ∈ en(m) is N P -firable from C if and only if it is firable (not necessarily N P -firable) from (m, IH(D)).

Proof. ⇐: trivial because IH(D) ⊆ D.

⇒: since t f is N P -firable from C, there exists an integer parameter valuation v, and ( θ, c, v) ∈ D such that for all transitions t i ∈ en(m), θ i ≥ θ f . And the result follows from Lemma 6 because v is an integer valuation.

Second, the cost of a class C = (m, D), for integer parameters, is cost

N (C) = inf ( θ,c,v)∈D,v∈N P c.
Lemma 8 is a direct consequence of Lemma 6: Proof. ⇒: if ( θ, c, v) ∈ IH(D) then it is also in D and the result follows from Lemma 1. ⇐: by Lemma 1, we know that there exists some ( θ, c, v) ∈ D, and since v is an integer valuation, by Lemma 6, ( θ, c, v) ∈ IH(D).

Lemma 10. Let C σ1 and C σ2 be two state classes such that IH(C σ1 ) IH(C σ2 ).

If a transition sequence σ is N P -firable from C σ1 it is also N P -firable from

C σ2 and cost N (C σ1.σ ) ≥ cost N (C σ2.σ ). Proof. Let C σ1 = (m 1 , D 1 ) and C σ2 = (m 2 , D 2 ). From Definition 3, for any point ( θ, c 1 , v) ∈ IH(D 1 ), there exists a point ( θ, c 2 , v) ∈ IH(D 2 ) such that c 2 ≤ c 1 .
With Lemma 7 and Lemma 8, this implies that: (i) cost N (C σ1 ) ≥ cost N (C σ2 ); (ii) if transition t is N P -firable from C σ1 , then it is N P -firable from C σ2 and Next(C σ1 , t) Next(C σ2 , t). And, as before, the result follows by a straightforward induction.

Algorithm 3 Restriction of (semi-)Algorithm 2 to integer parameter valuations. 1. if v ∈ PolyRes and v is an integer parameter valuation then there exists a run ρ in v(N ) such that cost(ρ) ≤ c max and lastm(ρ) ∈ Goal. 2. if Waiting = ∅ then for all integer parameter valuations v such that there exists a run ρ in v(N ) such that cost(ρ) ≤ c max and lastm(ρ) ∈ Goal, we have v ∈ PolyRes.

In Algorithm 3, we compute state classes as usual then handle them via their integer hulls. We can actually simply integrate integer hulls at the end of Algorithm 1 and use Algorithm 2 with this updated successor computation as proved by Lemma 11.

Termination of Algorithm 3

We now consider that parameter valuations are bounded by some value M 1 ∈ N (and that they still have integer values). We also assume that, for all integer parameter valuations, there exists M 2 ∈ Z such that for all runs ρ in v(N ), cost(ρ) ≥ M 2 : this allows us, as in [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF][START_REF] Bouyer | Symbolic optimal reachability in weighted timed automata[END_REF], to keep Algorithm 3 simple by doing away with negative cost loop-checking. Finally, we assume the net itself is bounded: there exists M 3 ∈ N such that for all reachable markings m, for all places p, m(p) ≤ M 3 .

To prove the termination of Algorithm 3 under these assumptions, we consider the symmetric relation to , such that x y iff y x. We prove that it is a well quasi-order (wqo), i.e., that for every infinite sequence of state classes, there exist C and C in the sequence, with C strictly preceding C such that C C . This implies that the exploration of children in Algorithm 3 will always eventually stop. Proposition 5. Let N be a bounded pcTPN, with bounded integer parameters and such that the cost of all runs is uniformly lower-bounded for all integer parameter valuations.

Relation is well-quasiorder on the set of state classes of N .

Proof. Consider an infinite sequence C 0 , C 1 , C 2 , . . . of state classes. Let C i = (m i , D i ).

From [START_REF] Boucheneb | Optimal reachability in cost time Petri nets[END_REF], we know that is a wqo for the state classes of bounded (non parametric) cost TPNs. So for each integer parameter valuation v, and using a classic property of wqo we can extract a subsequence of v(C 0 ), v(C 1 ), . . . that is completely ordered by . And since, we have a finite number of such parameter valuations, we can extract an infinite subsequence C i0 , C i1 , . . . such that for all integer parameter valuations v

, v(C i0 ) v(C i1 ) • • • .
Let us consider two of those: C ir and C is , with r < s.

Since IH(D is ) has integer vertices, and for any integer parameter valuation, v(C ir ) v(C is ), which implies that v(D is ) ⊆ v(D ir ), then all the vertices of D is are also in D ir . Now assume that some extremal ray r of D is is not in D ir . Then starting from some vertex x of D ir , there must be some λ ≤ 0 such that x + λr ∈ D is and the same holds for any λ ≥ λ (by convexity). But since r has rational coordinates for some value of λ , λ r is an integer vector and so is x + λ r, which contradicts the fact that v(D is ) ⊆ v(D ir ), for all integer parameter valuations v, and in particular ( x + λ r) |P . We can therefore conclude that D ir ⊆ D is and we now proceed to proving that D is is also "cheaper" than D ir .

We use another property of the vertices of convex polyhedra: vertices of a convex polyhedron of dimension n defined by m inequalities For the i-th lower-bound inequality on c, and the j-th lower-bound inequality on c , we define E ij as D in which we transform both constraints into equalities. Clearly, from the property above, this does not add any new vertex, but it may remove some. Second, by construction, we have

ij E ij = {( θ, min ( θ,c,v)∈IH(Di r ) c, min ( θ,c,v)∈IH(Di s ) c)| θ ∈ IH(D r ) |Θ }.
If we minimize c -c over E i , we know from the theory of linear programming that the minimum is obtained at a vertex of E ij , and therefore, in particular, for an integer valuation v of the parameters, and an integer vector θ of D ir . Since we have v(C ir ) v(C is ), we then know that for these values of the theta variables and parameters, c ≤ c . This means that this holds for the whole of E ij , and finally that C ir C is .

Case Study

We now consider a scheduling problem where some tasks include runnables, a key concept of the AUTomotive Open System ARchitecture (AUTOSAR), the open standard for designing the architecture of vehicle software [START_REF]Specification of RTE Software[END_REF]. Runnables represent the functional view of the system and are executed by the runtime of the software component [START_REF] Naumann | AUTOSAR runtime environment and virtual function bus[END_REF]. For their execution they are mapped to tasks and a given runnable can be split across different tasks to introduce parallelism, for instance. In industrial practice, runnables that interact a lot are mapped to the same task, in particular when they perform functions with the same period.

In this example, we consider 3 non-preemptive, periodic tasks T1, T2 and T3, on which have already been mapped some runnables that interact together; we add another independent runnable whose code must be split between tasks T1 and T2:

the period of task T1 is 100 time units; T1 includes a "fixed part", independent from the new runnable and whose execution lasts 22 t.u.; the period of T2 is 200 t.u.; T2 also has a fixed part lasting 28 t.u.; the period of T3 is 400 t.u.; its execution lasts 11 t.u.; the period of the runnable is 200 t.u.; its execution lasts 76 t.u.; parameter a denotes the duration of the section that is executed in T1 3 .

The processing unit consists of 2 cores C0 and C1; T3 can only execute on C0 whereas both T1 and T2 can execute on either core. When both cores are idle, the cost is null; when only one core is busy, the cost is equal to 2/t.u.; when both cores are busy, the cost is equal to 3/t.u. Any optimised strategy to divide the runnable over T1 and T2 and to allocate these tasks to C0 or C1 must therefore favour the cases where both cores are in the same state.

Figure 3 presents the model for this problem 4 . The associated cost function is: 2 * (C0 = C1) + 3 * C0 * C1 + 1000 * W 1 * (R1C0 + R1C1) + W 2 * (R2C0 + R2C1) + W 3 * R3C0 , where the name of a place (e.g. R1C0) represents its marking 5 .

We limit the study of the system to the first 400 t.u., at the end of which T1 has been executed 4 times, T2 twice and T3 once. A preliminary analysis (not detailed here for the sake of concision) showed that the lowest cost is 466. By setting our maximal cost to this value, we then check the following property with our Romeo tool: EF four==4 and two==2 and one==1 and cost≤466. The answer provided by Romeo is that the property is true iff a ∈ [START_REF] Hune | Linear parametric model checking of timed automata[END_REF][START_REF] Larsen | Model-Checking for Real-Time Systems[END_REF]. We then set a to 17; Romeo yields the following timed trace, in which the notation T1@t1 means that transition T1 is fired at date t1: T1C0@61, T2C1@69, T1@100, end1 C0@100, T1C0@100, end1 C0@139, end2 C1@139, T1@200, T2@200, T2C0@261, T1C1@261, T1@300, end1 C1@300, T1C1@303, end2 C0@331, T3C0@331, end3 C0@342, end1 C1@342 From this trace, we obtain the Gantt chart in Figure 4 (above). Setting a to 13 yields another timed trace, resulting in the Gantt chart in Figure 4 (below). In both cases, we can see that both cores are busy during 148 t.u. (and for 11 t.u., only one is idle), which confirms our analysis on the optimised strategy above.

3 Every 200 t.u., since T1 is executed twice as often as T2, T1 is running during (22 + a) * 2 = 44 + 2a t.u. whereas T2 is running during 28 + (76 -2a) = 104 -2a t.u. 4 To ensure a correct access to the cores, we could have added one place for each core and some arcs on each task to capture and release them but the resulting net would have been quite unreadable. Instead, we chose to add 2 integer variables C0 and C1 (both initialised to 0); a variable equal to 0 (resp. 1) obviously means the corresponding core is idle (resp. busy). 5 The last term ensures that such cases where an instance of a task is activated while a previous one is running are heavily penalised.

Conclusion

We have proposed a new Petri net-based formalism with parametric timing and cost features, thus merging two classic lines of work. For this formalism, we define an existential problem and a synthesis problem for parametric reachability within a bounded cost. We prove that the former is undecidable but we nonetheless give and prove a symbolic semi-algorithm for the latter. We finally propose a variant of the synthesis algorithm suitable for integer parameter valuations and prove its termination when those parameter valuations are bounded, and some other classic assumptions. This symbolic algorithm avoids the explicit enumeration of all possible parameter valuations. It is implemented in our tool Romeo and we have reported on a case-study addressing a scheduling problem, and inspired by the AUTOSAR standard. Further work includes computing the optimal cost as a function of parameters and investigating the case of costs (discrete and rates) as parameters.

  1.b) Initialise the parameters a and b such that: 0 < a ≤ 1 and 0 < b ≤ 1 1.a) Encoding of: when in state si, increment Cy and go to sj

Fig. 1 .

 1 Fig. 1. Increment Gadget (left) and Initial gadget (right)

1 .

 1 Either the machine halts, both counters C x and C y are bounded (let c their maximum value) and the halting and finite execution of the machine is within m steps. If c = 0 then the machine is a sequence of m zero-test taking m.b time units and the parametric Petri net N M can go within 1 time unit to a marking m halt if 0 < a ≤ 1 and 0 < b ≤ 1 m . If c > 0, since an instruction requires at most a + b time units, if a + b ≤ 1 m and if 0 < a ≤ b c then there exists a run that correctly simulates the machine, and eventually reaches m halt within 1 time unit. This set of valuations is non-empty: for example if c = 0, then we can choose a = b = 1 m and if c > 0, then, since m ≥ c, we can choose a = b m and b = 1 1+m hence a = 1 m(1+m) . 2. Or the machine does not halt. A step requires at least b time units then for any value v of the parameters, after a maximum number of steps (at most 1 b ), one whole time unit will elapse without v(N M ) reaching m halt .

  All the states of C σ share the same marking and can therefore be written as a pair (m, D) where m is the common marking

Lemma 1 .

 1 For all classes C σ = (m, D), ( θ, c, v) ∈ D if and only if there exists a run ρ in v(N ), and I : en(m) → I(Q ≥0 ), such that sequence(ρ) = σ, (m, I, c) = last(ρ), and θ ∈ I.

) 5 :

 5 eliminate (by projection) variables c, θi for all i from D 6: for all tj ∈ newen(m, t f ), add variable θ j to D , constrained by θ j ∈ Is(tj)Let C 0 = (m 0 , D 0 ) be the initial class. Domain D 0 is defined by the constraints ∀t i ∈ en(m 0 ), θ i ∈ I s (t i ), ∀t ∈ T, I s (t) = ∅, and c = 0. This gives a convex polyhedron of R |en(m0)|+|P|+1 ≥0

Definition 3 .

 3 Let C = (m, D) and C = (m , D ) be two parametric cost state classes. We say that C is subsumed by C , which we denote by C C iff m = m , D |Θ∪P ⊆ D |Θ∪P , and for all

Proof. 1 .

 1 By induction on the while loop: initially, PolyRes is empty so the result holds trivially. Suppose it holds after some iteration n, and consider iteration n + 1. Let v ∈ PolyRes after iteration n + 1. If v was already in PolyRes after iteration n then we can apply the induction hypothesis. Otherwise it means that if C σ = (m, D) is the class examined at iteration n + 1, then m ∈ Goal and v ∈ D ∩ (c ≤ c max ) |P . This means that there exists some point ( θ, c, v) ∈ D with c ≤ c max . By Lemma 1, this means that there exists a run ρ such that (m, I, c, v) = last(ρ), for some I such that θ ∈ I, and therefore lastm(ρ) ∈ Goal and cost(ρ) ≤ c max . 2. Let v be a parameter valuation such that there exists a run ρ in v(N ) such that cost(ρ) ≤ c max and lastm(ρ) ∈ Goal. Let σ = sequence(ρ). We proceed by induction on the length n of the biggest suffix σ 2 of σ such that, either σ 2 is empty or, if we note σ = σ 1 σ 2 , with the first element of σ 2 being transition t, then C σ1t ∈ Passed. If n = 0, then C σ = (m, D) ∈ Passed. By Lemma 1, v ∈ D |P and m ∈ Goal. From the latter, we have D ∩ (c ≤ c max ) |P ⊆ PolyRes and therefore v ∈ PolyRes because v ∈ D ∩ (c ≤ c max ) |P . Consider now n > 0 and assume the property holds for n -1. Since n > 0, then there exists a transition t and a sequence σ 3 such that σ 2 = t.σ 3 . By definition of σ 2 , we have C σ1 ∈ Passed but C σ1.t ∈ Passed. By Lemma 4, since Waiting = ∅, there must exists some class C σ such that C σ1.t C σ .

Lemma 6 .

 6 Let (m, D) be a state class of a pcTPN and let ( θ, c, v) be a point in D. If v is an integer valuation, then ( θ, c, v) ∈ IH(D). Proof. Since ( θ, c, v) ∈ D then ( θ, c) ∈ v(D). By Proposition 3, v(D) being the firing domain of a state class in a (non-parametric) cost TPN, it has integer vertices, and therefore v(D) = IH(v(D)). Point ( θ, c) is therefore a convex combination of integer points in v(D). Clearly, for all integer points ( θ , c ) in v(D), we have that ( θ , c , v) is an integer point of D. Since D is convex, this implies that ( θ, c, v) ∈ IH(D).

Lemma 8 .Lemma 9 .

 89 Let (m, D) be a state class. We have: cost N ((m, D)) = cost((m, IH(D)). If v is an integer parameter valuation, then for all classes C σ = (m, D), ( θ, c, v) ∈ IH(D) if and only if there exists a run ρ in v(N ), and I : en(m) → I(Q ≥0 ), such that sequence(ρ) = σ, (m, I, c) = last(ρ), and θ ∈ I.

1: PolyRes ← ∅ 2 :Proposition 4 .

 24 Passed ← ∅ 3: Waiting ← {(m0, D0)} 4: while Waiting = ∅ do 5: select Cσ = (m, D) from Waiting 6: if m ∈ Goal then 7: PolyRes ← PolyRes ∪ IH(D) ∩ (c ≤ cmax) |P 8: end if 9: if for all C ∈ Passed, IH(Cσ) IH(C ) then 10: add Cσ to Passed 11: for all t ∈ firable(IH(Cσ)), add Cσ.t to Waiting 12: end if 13: end while 14: return PolyResUsing Lemma 9 instead of Lemma 1, and Lemma 10 instead of Lemma 3 in the proof of Proposition 2, we get the following proposition, stating the completeness and soundness of Algorithm 3. After any iteration of the while loop in Algorithm 3:

Lemma 11 .

 11 Let (m, D) be a state class of a pcTPN N , and t a transition firable from C. Let (m , D ) = Next((m, D), t) and (m , D ) = Next((m, IH(D)), t). Then m = m and IH(D ) = IH(D ). Proof. The equality of markings is trivial so we focus on firing domains. By definition of the integer hull, we have IH(D) ⊆ D. Since the computation of the next class domain is non-decreasing with respect to inclusion, we then have D ⊆ D . Taking the integer hull is also non-decreasing wrt. inclusion, so IH(D ) ⊆ IH(D ). Consider now an integer point ( θ , c , v) in D . Then ( θ , c ) ∈ v(D ). Consider state class computations in the (non-parametric) cost TPN v(N ): there exists some point ( θ, c) in v(D) such that (m , θ , c ) ∈ Next((m, {( θ, c)}), t). Since ( θ, c, v) thus belongs to D and since v is an integer parameter valuation, by Lemma 6, we have that ( θ, c, v) ∈ IH(D). Thus ( θ , c , v) ∈ D and since it is an integer point, it is in IH(D ).

  n k=1 a kl x k ≤ b l , for j ∈ [1..m] are solutions of a system of n linearly independent equations n k=1 a kl x k = b l , with l in a subset of size n of [1..m]. Now consider the polyhedron D obtained from IH(D ir ), with its cost variable c, by adding one variable c constrained by the cost inequalities of IH(D is ). Clearly, since c and c are not constrained together, the vertices of D are those of IH(D ir ), extended with the corresponding minimal and maximal values of c , and symmetrically those of IH(D is ), extended with the corresponding minimal and maximal values of c . Since the inequalities constraining c and c have integer coefficients, and IH(D is ) and IH(D ir ) have integer vertices, D also has integer vertices.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Offline non preemptive scheduling problem

  line 11) and C σ ∈ Passed (line 10). Since we add at most one state class to Passed at each iteration, C σ was added in a previous iteration and we can apply to it the induction hypothesis, which allows us to prove the first part of the invariant; 2. lines 6 and 7 obviously imply the second part of the invariant; 3. if C σ ∈ Passed, then the condition of the if on line 9 is true and then for any transition t that is firable from C σ , C σ.t is added to Waiting (line 11) so the third part of the invariant holds for C σ . Nevertheless, C σ itself is no longer in Waiting, and it is (except for the initial state class) the successor of some state class in Passed. But then we have only two possibilities: either C σ has been added to Passed in line 10 if the condition on line 9 was true, and certainly C σ C σ , or there exists C ∈ Passed such that C σ C if that condition was false. Therefore the third part of the invariant holds.

This work is partially supported by the ANR national research program PACS (ANR-14-CE28-0002)


J A A Lowther
email: james.lowther@cefas.co.uk

A Bosch

S Butot

J Ollivier

D Mäde

S A Rutjes

G Hardouin

B Lombard

P In't Veld

Alexandre Leclercq

Validation of EN ISO method 15216 -Part 1 -Quantification of hepatitis A virus and norovirus in food matrices

Keywords: norovirus, hepatitis A virus, real-time RT-PCR, validation, standardisation, bivalve molluscan shellfish, soft fruit, vegetables, bottled water, surfaces

Introduction

Hepatitis A virus (HAV) and norovirus are important agents of food-borne human viral illness. The foodstuffs most commonly linked to illnesses caused by these viruses include bivalve molluscan shellfish (reviewed in [START_REF] Bellou | Shellfish-borne viral outbreaks: a systematic review[END_REF], soft fruit [START_REF] Bernard | Large multistate outbreak of norovirus gastroenteritis associated with frozen strawberries, Germany[END_REF][START_REF] Severi | Europe associated with consumption of frozen berries, 2013 to 2014[END_REF] and a variety of leaf, stem or bulb vegetables including lettuce [START_REF] Ethelberg | Outbreaks of gastroenteritis linked to lettuce, Denmark[END_REF] and green onions [START_REF] Dentinger | An outbreak of hepatitis A associated with green onions[END_REF]. Separately the presence of norovirus RNA in bottled water has been variously reported and debated [START_REF] Beuret | Norwalk-like virus sequences detected by reverse transcriptionpolymerase chain reaction in mineral waters imported into or bottled in Switzerland[END_REF][START_REF] Blanco | Norovirus in bottled water associated with gastroenteritis outbreak, Spain[END_REF][START_REF] Sanchez | Presence of norovirus sequences in bottled waters is questionable[END_REF]. Outbreaks of viral illness due to contamination of the surfaces of foods, or food preparation surfaces by for example infected food handlers have also been documented [START_REF] Chen | An outbreak of norovirus gastroenteritis associated with asymptomatic food handlers in Kinmen, Taiwan[END_REF][START_REF] Thornley | Multiple outbreaks of a novel norovirus GII.4 linked to an infected post-symptomatic food handler[END_REF]. Virus analysis can contribute towards risk management for these hazards. However, until recently no standard method has existed for virus analysis in foods and it is documented that different methods can give divergent results (Lees and CEN WG6 TAG4, 2010). For these reasons a European project to develop a standardised method to detect these viruses in a variety of food matrices was launched in 2004 by CEN/TC275/WG6, the European Committee for Standardisation (CEN) Working Group on the Microbiology of the Food Chain, with a Technical Advisory Group convened specifically for this purpose (CEN/TC275/WG6/TAG4). Method development within this group proceeded by a combination of consensus, methodological ring-trials within the group and methodological investigations by individual group members. This project culminated in the publication in 2013 of a two-part technical specification for determination of the viruses in food matrices in 2013 (ISO/TS 15216:2013;Anonymous, 2013a and[START_REF] Anonymous | Microbiology of food and animal feed --Horizontal method for determination of hepatitis A virus and norovirus in food using real-time RT-PCR --Part 2: Method for qualitative detection[END_REF].

As part of the European Mandate No. M381 to validate 15 standards in the field of food microbiology, this study aimed to validate ISO 15216 -Part 1, the method for quantification, in seven food matrices; bottled water, food surfaces (bell pepper pieces), Pacific oysters (Crassostrea gigas), common mussels (Mytilus edulis), raspberries, lettuce and green onions, in order to replace the technical specification with a full, validated EN/ISO standard.

Materials and methods

The method evaluated

As no routine methods exist for culture of norovirus, and HAV culture methods [START_REF] Flehmig | Hepatitis A-virus in cell culture: I. propagation of different hepatitis A-virus isolates in a fetal rhesus monkey kidney cell line (Frhk-4)[END_REF] are not appropriate for routine application to food matrices, detection is reliant on molecular methods using the reverse-transcriptase polymerase chain reaction (RT-PCR). In this study, testing of all samples for norovirus GI, norovirus GII and HAV followed matrix-specific test protocols compliant with ISO 15216-1, using methods for RNA extraction and real-time RT-PCR, real-time RT-PCR primers and probes [START_REF] Costafreda | Development, evaluation, and standardization of a realtime TaqMan reverse transcription-PCR assay for quantification of hepatitis A virus in clinical and shellfish samples[END_REF][START_REF] Da Silva | Evaluation of removal of noroviruses during wastewater treatment, using real-time reverse transcription-PCR: different behaviors of genogroups I and II[END_REF][START_REF] Hoehne | Detection of Norovirus genogroup I and II by multiplex real-time RT-PCR using a 3'-minor groove binder-DNA probe[END_REF][START_REF] Loisy | Real-time RT-PCR for norovirus screening in shellfish[END_REF][START_REF] Svraka | Etiological role of viruses in outbreaks of acute gastroenteritis in The Netherlands from 1994 through 2005[END_REF] and the process control virus (mengo virus strain MC0; [START_REF] Costafreda | Development, evaluation, and standardization of a realtime TaqMan reverse transcription-PCR assay for quantification of hepatitis A virus in clinical and shellfish samples[END_REF], detailed in the informative annexes. For each matrix, the method consisted of a matrix-specific virus extraction followed by common RNA extraction and real-time RT-PCR detection elements. Briefly, for food surfaces, virus extraction used swabbing of the surface with a sterile cotton swab, followed by elution into lysis buffer. For fruit and vegetable matrices, virus extraction was by elution with agitation followed by precipitation with PEG/NaCl (with additional extraction steps for pectin rich fruit). For bottled water, adsorption and elution using positively charged membranes followed by concentration by ultrafiltration was used and for bivalve molluscan shellfish, viruses were extracted from the tissues of the digestive glands using treatment with a proteinase K solution. Virus extracts from all matrices were subjected to a common RNA extraction method based on virus capsid disruption with chaotropic reagents followed by adsorption of RNA to silica particles. Detection of virus sequences within the sample RNA utilised real-time RT-PCR with hydrolysis probes in duplicate reactions for each sample RNA and target virus combination. Undiluted and 1/10 diluted RNA was tested for each sample; in accordance with ISO 15216-1, the results for undiluted RNA were used unless this demonstrated significant RT-PCR inhibition, in which case results for 1/10 diluted RNA would be checked. Due to the complexity of the method, a comprehensive suite of controls was included, including negative process, extraction and RT-PCR controls, and controls for RT-PCR inhibition and extraction efficiency. Quantification of target copies per microliter of sample RNA was by reference to a standard curve generated from a dilution series of dsDNA carrying the relevant target sequence as described in Annex G of ISO 15216-1:2017 (Anonymous, 2017b). In accordance with ISO 15216-1, quantities were not corrected according to extraction efficiency and RT-PCR inhibition results.

Design of the study

For each matrix, the study comprised two parts; part 1; method characterisation and part 2; interlaboratory study. For each matrix part 1 was carried out by a single expert laboratory (Table 1) testing 60 samples of matrix contaminated with a dilution series of the three target viruses in order to determine method characteristics including limit of detection and limit of quantification. Part 2 was carried out by one organising expert laboratory and 10 participating laboratories per matrix. Each participating laboratory tested duplicate samples, prepared and distributed by the organising laboratory, designed to represent 4 levels of contamination (high, medium, low, and negative) in order to determine method characteristics including reproducibility and repeatability. In total, 18 laboratories from 11 countries in Europe participated in one or both parts of the study for one or more matrices.

In its role as project leader, the European Union Reference Laboratory for Monitoring Bacteriological and Viral Contamination of Bivalve Molluscs, based at the Centre for Environment, Fisheries and Aquaculture Science, Weymouth, United Kingdom, was responsible for management of the validation project, development of study protocols, generation and distribution of control materials and virus stocks, and collation and analysis of study data.

Part 1; method characterisationgeneral considerations

For part 1 of the study, test sample sizes were 2g (bivalve shellfish digestive tissues), 25g (fruit and vegetables), 330ml (bottled water) and 25cm 2 (bell pepper pieces to represent the food surfaces matrix).

All samples were tested using the relevant matrix-specific test protocols (compliant with ISO 15216-1) then both calculated quantities and raw data were forwarded to the project leader for quality checking according to ISO 15216-1 and analysis as described under Generation of Method Characteristics below.

Part 1; method characterisationpreparation of the virus dilution series

The viruses used for contamination were genotypes GI.4 and GII.4 norovirus from faecal suspensions and HM175/43c strain HAV derived from tissue culture. The limited availability of norovirus stocks (particularly GI) meant that a simple log10 dilution series with sufficiently high starting levels to allow for contamination across a wide range of dilutions was not practical; for this reason a 0.5 log10 dilution series (~3.16x) was used instead. For each matrix 0.5 log10 dilution series of contaminated matrix samples at 9 separate levels were prepared as below (different strategies reflect the modes of natural contamination of the different food matrices). In each case virus dosing was calculated to provide detectable levels of 100 -1,000 copies/μl in the RNA extract at the highest contamination level (equivalent to ~20,000 -200,000 copies/g for bivalve shellfish samples, 400 -4,000 copies/g for fruit and vegetable samples, 400 -4,000 copies/cm 2 for food surface samples and 30 -300 copies/ml for bottled water samples), based on the results of a trial contamination:a)

Bivalve shellfish matrices; contaminated matrix was prepared by bioaccumulation of shellfish with the 3 target viruses. Food containing the viruses was added to a tank containing shellfish undergoing filter feeding behaviour, then shellfish were harvested after 18hrs

(oysters) or 24hrs (mussels). Contaminated digestive tissues were dissected and homogenised by blending, then a ≈0.5 log10 dilution series of this contaminated matrix in homogenised uncontaminated matrix (digestive tissues from clean shellfish) was prepared by serial blending of 3.8 g of contaminated matrix in 8.2 g of uncontaminated matrix.

b) Non-shellfish matrices; a mixture of the 3 target viruses at high initial levels was prepared in buffer, then a ≈0.5 log10 dilution series of this was prepared by serial dilution of e.g. 100 μl contaminated virus mix in 216 μl buffer. Clean matrix samples were then contaminated with the different dilutions. For fruit, vegetable and food surface samples, the virus solution was pipetted across the surface of the samples, then left to air dry for 20 minutes in a laminar flow cabinet. For bottled water samples, the virus solution was added directly to the sample.

Part 1; method characterisationstructure of the study

For each matrix a total of 54 contaminated samples were prepared as follows (the different strategies reflect the time practicalities of virus extraction methods for the different food matrices; strategies a) and b) with triplicate subsamples were used where practicable, strategies c) or d) were used where it was not practical to perform virus extraction on 27 samples simultaneously):a) Bivalve shellfish matrices; triplicate subsamples were prepared for each level of the 0.5 log10 dilution series as described above. This entire procedure was carried out on two separate occasions (9 dilution levels x 3 subsamples per level x 2 occasions = 54 samples).

b) Food surface matrix; for each level of the 0.5 log10 dilution series of virus mix triplicate subsamples of clean matrix (bell pepper pieces) were prepared. In each case, the virus solution was pipetted across the exterior surface of the samples, then left to air dry for 20 minutes in a laminar flow cabinet. This entire procedure was carried out on two separate occasions (9 dilution levels x 3 subsamples per level x 2 occasions = 54 samples). c) Bottled water matrix; single samples were prepared for each level of the 0.5 log10 dilution series of virus mix. In each case, the virus solution was added directly to the sample.

This entire procedure was carried out on six separate occasions (9 dilution levels x 6 occasions = 54 samples).

d)

Fruit and vegetable matrices; single samples were prepared for each level of the 0.5 log10 dilution series of virus mix. In each case, the virus solution was pipetted across the exterior surface of the samples, then left to air dry for 20 minutes in a laminar flow cabinet.

This entire procedure was carried out on six separate occasions (9 dilution levels x 6 occasions = 54 samples).

In addition to 54 contaminated samples, 6 samples of uncontaminated matrix were tested in parallel to give a total of 60 test samples per matrix.

Part 2; interlaboratory studygeneral considerations

For part 2 of the study, test sample sizes were 2g (bivalve shellfish digestive tissues), 25g (fruit and vegetables), 300ml (bottled water) and 50cm 2 (bell pepper pieces to represent the food surfaces matrix).

For each matrix interlaboratory study all samples were tested using the relevant matrix-specific test protocols (compliant with ISO 15216-1) then both calculated quantities and raw data were forwarded to the project leader for quality checking according to ISO 15216-1 and analysis as described under Generation of Method Characteristics below.

Part 2; interlaboratory studypreparation of test samples

The viruses used for contamination were the same as for part 1. For each of the 7 matrices the relevant expert laboratory prepared four batches of test samples each comprising multiple samples containing high, medium, low or negative levels of each of the three target viruses. In each case virus dosing was calculated to provide detectable levels of 20 -200 copies/μl in the RNA extract at the highest contamination level (equivalent to ~4,000 -40,000 copies/g for bivalve shellfish samples, 80 -800 copies/g for fruit and vegetable samples, 40 -400 copies/cm for food surface samples and 6.7 -67 copies/ml for bottled water samples), based on the results of a trial contamination. For medium and low levels respectively, intended levels were 1/5 th and 1/25 th those in the high levels samples respectively.

For Pacific oysters a single bioaccumulation to produce highly contaminated digestive tissues was carried out. This material was homogenised then high, medium and low level batches were prepared by diluting this contaminated matrix in uncontaminated matrix (digestive tissues from uncontaminated Pacific oysters) to the appropriate levels. Through this approach it was aimed to produce homogenous starting materials with a clearly defined proportional difference between the different contamination levels. For common mussels due to the practical difficulties of dissection of sufficient quantities of digestive tissues in a single laboratory (the digestive tissues of individual mussels are smaller than for oysters), three separate bioaccumulations to produce high, medium and low levels were carried out.

For other matrices, multiple portions of matrix were directly contaminated with virus stocks containing a mix of all three target viruses at high, medium or low levels using contamination methods as described for part 1.

To demonstrate adequate homogeneity and stability a minimum of ten samples were tested for each contamination level by the expert laboratory using the relevant matrix-specific protocol. The number of samples and testing schedule was dependent on the storage temperature and the required lifetime of test samples as described below.

Part 2; interlaboratory studytesting by participant laboratories

For each matrix eight anonymised test samples (duplicate samples for each contamination level) were sent to each of ten participating laboratories (expert laboratories did not participate in the interlaboratory study for matrices where they prepared test samples). Each participating laboratory tested the eight samples using the relevant matrix-specific protocol within a specified timescale. Depending on whether it was practical to freeze samples (due to considerations including damage to the samples from freezethaw cycles) test samples were either distributed frozen or chilled. For matrices where test samples were distributed frozen (bivalve molluscs, green onions, raspberries) laboratories were permitted to store samples frozen prior to testing, but were instructed to return data to the project leader by a specified date 2 months from the receipt of samples. For matrices where test samples were distributed chilled (lettuce, bottled water, food surfaces) laboratories were instructed to carry out virus and RNA extractions within 48hrs of receipt, and were further instructed to return data to the project leader by a specified date 1 month from the receipt of samples. In all cases protocols required participating laboratories to respect maximum storage times and temperatures of intermediate test materials (virus extract and RNA).

Generation of method characteristicsgeneral considerations

Results for test samples that were invalid according to the protocols due to unacceptable extraction efficiency or RT-PCR inhibition levels were excluded from the analysis, as were false negative results.

Obtained results below the theoretical limit of detection (tLOD; the concentration equivalent to the detection of a single target copy across the two real-time RT-PCR reactions) were adjusted upwards to the tLOD.

Generation of method characteristicslimit of detection and limit of quantification

Limit of detection (LOD) and limit of quantification (LOQ) characteristics were determined using the data generated in part 1; method characterisation. For these data sets, "anticipated results" (the designated correct result) for each contaminated test sample were calculated (separately for each dilution series)

as follows:a)

Bivalve shellfish matrices and food surfaces (matrices where triplicate subsamples at each level of the dilution series were used); the anticipated result at each level for each dilution series was calculated as the geometric mean of the obtained results for the 3 subsamples at the highest concentration within the series (≡ anticipated result at the highest concentration) multiplied by the dilution factor. For example, for oysters, norovirus GI, occasion 1, obtained results for the 3 subsamples at the highest concentration were 18036 copies/g, 21220 copies/g and 9415 copies/g respectively (geometric mean = 15331 copies/g). The anticipated results at the different dilutions of the series were therefore calculated at 15331, 4848 ( 15331x 10 -0.5 ), 1533, 484, 153, 48, 15, 5 and 2 copies/g respectively. b) Bottled water, fruit and vegetable matrices (matrices where single samples at each level of the dilution series were used); the anticipated result at each level for each dilution series was calculated as the antilog of the intercept of the line of best fit for the plot of log10 positive obtained results vs. log10 dilution factor (≡ anticipated result at the highest concentration within the series), multiplied by the dilution factor. For example, for bottled water, norovirus GI, occasion 1, results for the 9 samples within the dilution series were 141.84 copies/ml, 56.86 copies/ml, 19.72 copies/ml, 4.35 copies/ml, 0.75 copies/ml, 0.24 copies/ml, 0.08 copies/ml, 0.04 copies/ml and not detected. The intercept of the line of best fit of the log10 positive obtained results vs. log10 dilution factor was 2.23, corresponding to an antilog of 170.01 copies/ml. The anticipated results at the different dilutions of the series were therefore calculated at 170.01, 53.76 (170.01 x 10 -0.5 ), 17.00, 5.37, 1.70, 0.54, 0.17, 0.05 and 0.02 copies/ml respectively.

The LOD and LOQ were then calculated for each matrix/target virus combination using the datasets of anticipated vs. obtained results as follows:a) the LOD95 (the lowest concentration of target virus that can be consistently detected in 95% of samples tested under routine laboratory conditions) was calculated using an online

Microsoft EXCEL program for the estimation of the POD (probability of detection) function and the LOD of a qualitative microbiological measurement method according to [START_REF] Wilrich | Estimation of the POD function and the LOD of a qualitative microbiological measurement method[END_REF] using the number of positive and negative results at each anticipated level (http://www.wiwiss.fu-berlin.de/fachbereich/vwl/iso/ehemalige/wilrich/PODLOD_ver9.xls accessed 10th November 2017).

b) the LOQ was calculated for each matrix/target virus combination using a method adapted from [START_REF] Armbruster | Limit of Blank, Limit of Detection and Limit of Quantitation[END_REF]. This characteristic was determined by looking at the log10 transformed data from the half log10 immediately below the LOD and higher (e.g. if LOD was 1.85 log10 then results from an anticipated level of 1.5 log10 and upwards were used). A regression line was fitted to these selected data and the residuals calculated (differences between observed value and fitted line). The standard deviation of these residuals was calculated in one log10 intervals of the data (moving up in half log10 steps) and the LOQ was determined as the level above which the standard deviation was always below 0.5 log10. The LOQ was set no lower than the LOD.

Generation of method characteristicsrepeatability and reproducibility

Repeatability and reproducibility characteristics were determined using the data generated in part 2;

interlaboratory study. All data returned by the participant laboratories was log10 transformed then outlying results for each matrix/target virus/contamination level combination were identified using Mandel's h and k statistic (test for graphical consistency; [START_REF] Mandel | A new analysis of interlaboratory test results[END_REF], Cochran's test for withinlaboratory variability [START_REF] Cochran | The distribution of the largest of a set of estimated variances as a fraction of their total[END_REF] and Grubbs' test for between-laboratory variability [START_REF] Grubbs | Sample criteria for testing outlying observations[END_REF]. Pairs of results identified by these tests as outliers were removed from the data set provided this did not result in fewer than 16 valid, positive data points remaining in the set (out of 20 total). Where individual laboratories recorded a high proportion (>50%) of invalid, false negative or outlying results across a given matrix, the entire data set for that laboratory/matrix combination (for all target viruses and all contamination levels) was excluded from the analysis. Repeatability and between-laboratory variances were calculated according to the formulae provided in ISO 5725-2:1994 (Anonymous, 1994).

The reproducibility variance was then calculated as the sum of the repeatability and between-laboratory variances, and the repeatability and reproducibility standard deviations determined as the square root of the respective variances. Repeatability and reproducibility limits were calculated as the respective standard deviations multiplied by a fixed factor of 2.8 (~2√2). These limits are defined as the absolute difference between two independent single (log10-transformed) test results or the ratio of the higher to the lower of the two test results on the normal scale, obtained under repeatability or reproducibility conditions respectively, that will not be exceeded in more than 5% of cases.

Results and Discussion

Recalibration of norovirus GII results

After completion of all practical work it became apparent through unconnected sequencing analysis carried out by a CEN/TC275/WG6/TAG4 group member that the dsDNA quantification standard for norovirus GII used in both parts of the validation comprised a mix of two distinct sequences (sequence y and sequence zsee Figure 1), rather than the intended sequence (sequence x). Sequences y and z each included mismatches relative to the primer and probe set used (major for sequence y, including several missing nucleotides in the probe target region, minor for sequence z). As a result, the mixed dsDNA standard exhibited reduced amplification efficiency and increased Cq values. This caused an upward bias in quantitative results obtained for the test samples relative to those that would have been obtained using the intended sequence.

To rectify this issue a calibration factor was used to adjust GII test sample results to compensate for the quantification bias introduced by the incorrect sequences. To establish the appropriate calibration factor, the GII dsDNA quantification standard used in the validation (adjusted to 1x10 5 copies/µl following quantification by spectrophotometry) was subjected to qPCR alongside a dilution series prepared using newly synthesised dsDNA quantification control with no mismatches relative to the primer and probe sequences (sequence x). This established that the validation standards exhibited an effective 3.157x consistent reduction in amplification efficiency (≈1.66 Cq values), resulting in a corresponding 3.157x overestimation in quantitative results for test samples in the validation (results not shown). All GII results in both parts of the validation were therefore adjusted downwards by this same factor. This issue affected absolute quantification but not relative quantification of samples within and between laboratories. Therefore, of the performance characteristics determined, only LOD and LOQ were affected by this adjustment. Repeatability and reproducibility characteristics were not impacted.

Specific investigations to identify the root cause of the mixed sequences in the GII dsDNA quantification standard were not carried out. However, it seems likely that sequence errors were introduced either during the synthesis of the insert sequence, or the cloning of the insert into the plasmid vector. The presence of two sequences in the purified plasmid preparation could be explained if during the cloning procedure a mixed colony of transformed E. coli host cells containing two different clonal populations (due to two different transformant cells colonising the selection media in close proximity) was inoculated into the growth medium used for plasmid preparation.

As a result of this issue a mandatory sequence verification for quantification standards (either by the user laboratory or the supplier) was introduced to the text of ISO 15216-1:2017 (Anonymous, 2017b) to ensure quantification standards contain the intended sequences.

Dilution of quantification standards using water only

Quantification standards prepared for part 1 of the validation study were diluted to working concentration, and further diluted to produce the standard curve, using molecular biology grade water as the diluent, consistent with the text of ISO/TS 15216-1:2013 (Anonymous, 2013a). During the course of part 1 it became apparent that sporadic problems with stability and homogeneity of standards, and the production of the standard curve, had occurred. Analysis of quantification standard curve data from part 1 of the validation study (data not shown) indicated that significant issues with the performance of the quantification standards were experienced in the studies on the food surfaces and green onions matrices, with problems including high intercept values, significant variability in intercept values between real-time RT-PCR runs and poor PCR efficiency values (slope <-3.6). For other matrices no such issues were apparent. For these reasons LOD and LOQ values generated for the affected matrices in part 1 were considered unreliable and unrepresentative of the affected matrices and are not included here.

Prior to the start of part 2 of the validation a series of experiments to investigate this issue were carried out (data not shown). It was determined that the issues experienced had been caused by the use of water as diluent in the dilution of quantification standards to working concentration and for the production of the standard curve, and that these issues could be eliminated by using an appropriate buffer (e.g. TE buffer) as diluent. Accordingly, this practise was adopted for part 2 of the validation study and no problematic results were noted. As a result of these findings, the text of ISO 15216-1:2017 (Anonymous, 2017b) specifically mandates the use of an appropriate buffer for such dilutions, rather than water only.

RT-PCR inhibition in the common mussel matrix

During part 1 of the validation study for the common mussel matrix a very high proportion (64.4%) of test samples provided results above the acceptable RT-PCR inhibition threshold (75% as detailed in ISO 15216-1) when undiluted RNA was tested. A smaller proportion (11.1%) also provided results above the acceptable threshold when 1/10 diluted RNA was tested. Accordingly, the majority of test sample results for the mussel matrix were determined using diluted RNA, as per the text of ISO 15216-1. The protocol for bivalve mollusc matrices for part 1 of the validation study specified the preparation of a dilution series from a single stock of contaminated digestive tissues diluted in a single stock of uncontaminated digestive tissues. It was not possible within the study protocol to take account of possible sample-to-sample matrix variations affecting parameters such as RT-PCR inhibition. From this study it is not clear whether high levels of RT-PCR inhibition are commonly found in common mussel samples, or whether the issue was specific to the particular sample of mussels selected for use in the study. Of note however in part 2 of the validation study only 1.8% of results for common mussel samples were affected by high RT-PCR inhibition levels. It can reasonably be assumed that the use of results from diluted RNA for the majority of test samples in part 1 will have impacted the LOD and LOQ values calculated for common mussels, and that these method characteristics may be unrepresentative of the matrix as a result. For this reason, they are not included here. Further data would therefore be required using non-inhibitory mussels to establish representative LOD and LOQ values for this matrix.

Results of part 1; method characterisation

As an example of the data obtained in part 1; method characterisation, the results for norovirus GI in the bottled water are shown as a plot of obtained versus anticipated results in Figure 2. The probability of detection/limit of detection plot derived from this data according to [START_REF] Wilrich | Estimation of the POD function and the LOD of a qualitative microbiological measurement method[END_REF] is shown in Figure 3. For this dataset the determined LOD was 0.18 copies/ml. For determination of LOQ a lineof-best-fit was prepared for all data points corresponding to log10 transformed anticipated levels above 0.1 copies/ml (shown in Figure 2). The standard deviations of the residuals from the line-of-best fit of all data points in one log10 intervals are shown in Table 2. In all cases these were below 0.5; the LOQ for norovirus GI in bottled water was therefore set at the same level as the LOD. Results were generated for the other study matrices using the same approach (calculations not shown).

Results of part 2; interlaboratory study

As an example of the data obtained in part 2; interlaboratory study, the results for norovirus GI in the bottled water for the high, medium, low and negative levels are shown in Figure 4. Test samples for this matrix were distributed chilled; the expert laboratory tested 5 subsamples at each level prior to the distribution ("post-preparation") and five subsamples at the end of the period for testing allotted to the participating laboratories ("post distribution"). Following quality control checking, one negative result (Lab 08, low, subsample 1) was removed from the data set. The method characteristics derived from the dataset for this matrix/target virus combination are shown in Table 3. Results were generated for the other study matrices using the same general approach (calculations not shown).

Method characteristics for all matrices

The LOD, LOQ, repeatability standard deviations and reproducibility standard deviations calculated for the seven matrices under examination are given in Tables 4 to 10. The repeatability and reproducibility standard deviations shown are the averages of those obtained at high, medium and low contamination levels; separate repeatability and reproducibility standard deviations and limits for all contamination levels are provided in Annex J of ISO 15216-1:2017 (Anonymous, 2017b).

Conclusion

This paper details the validation study carried out on ISO 15216-1, the standard method for quantification of hepatitis A virus and norovirus in foods. Detailed information on the performance of the method in seven different matrices was generated and has been included in ISO 15216-1:2017 (Anonymous, 2017b), the newly published edition of the standard. Acceptability criteria for the method characteristics have not been determined. The values obtained for repeatability and reproducibility in this study are comparable to other quantitative (enumeration) methods for bacterial pathogens validated at the same time under European Mandate No. M38 however. For example, ISO 10272-2, the colony count method for Campylobacter spp. (Anonymous, 2017a) produced average repeatability and reproducibility standard deviations of 0.20 and 0.40 log10 respectively, compared with 0.23 and 0.50 for the virus method. For LOD and LOQ characteristics it was notable that in the majority of cases the values determined for HAV were higher than for norovirus GI and GII. The primer set used in this study amplifies a relatively long product (157-188 bp depending on the strain of HAV), while for norovirus shorter, more optimal products are amplified (86 bp and 89 bp for GI and GII respectively; Anonymous, 2017b). This difference may have accounted for the different relative sensitivities of the methods for the different target viruses. Design of a broadly-reactive real-time RT-PCR primer probe set for HAV that amplifies a shorter product is complicated by the sequence diversity across strains however.

A number of technical issues were encountered through the validation study. It was necessary to recalibrate results for norovirus GII due to an issue with the sequence of the quantification standard. A mandatory sequence check has been added to the text of the new ISO, however the potential issue of variability between plasmids prepared in different laboratories remains. In this context, the availability of verified standard materials would be desirable. The European Union Reference Laboratory for Monitoring Bacteriological and Viral Contamination of Bivalve Molluscs supports laboratories within its network through the provision of quantification standards, however introduction of commercially available standards would be beneficial, to reduce one possible source of variability in results between laboratories.

For the green onions and food surfaces matrices it was not possible to determine LOD and LOQ characteristics due to problems with the performance of the quantification standards. These problems were completely alleviated through the use of TE buffer to dilute the standards; on this matter the literature is diverse with many contradictory examples of the use of either water or a variety of different buffers as a diluent, and no clear recommendations within journal publications in either direction. This presumably reflects the sporadic and unpredictable nature of the problems encountered, as in many cases the use of water will not result in problems (and indeed did not always in this validation study).

However, the introduction of a mandatory requirement for dilution of standards with a suitable buffer should ensure more stable performance.

Finally, problems were encountered during part 1 of the validation for mussels, due to the inhibitory nature of the matrix. This raises questions of the universal applicability of the method described in the ISO to different types of foods within the broad categories of bivalve molluscs, leaf, stem and bulb vegetables etc. In the case of mussels it seems likely that the problems encountered were due to an unusually inhibitory batch; many TAG4 members have successfully applied the ISO method to the analysis of mussels. However, the possibility of certain species of bivalve or types of vegetable providing consistently poor results remains; laboratories encountering such issues should consider whether the method in ISO-15216 is appropriate in these cases.

In addition to ISO 15216-1, a second part to the standard, ISO 15216-2, detailing a method for qualitative detection (not quantification) was published as a technical specification in 2013 [START_REF] Anonymous | Microbiology of food and animal feed --Horizontal method for determination of hepatitis A virus and norovirus in food using real-time RT-PCR --Part 2: Method for qualitative detection[END_REF]. Although the validation study was not designed to examine this part of the method, the data generated is sufficient to determine relevant method characteristics such as LOD50, specificity and sensitivity, and CEN/TC275/WG6 and its ISO sister group, ISO/TC34/SC9 have resolved that a new revision of ISO 15216-2, harmonised with ISO 15216-1:2017, and including method characteristics determined in this way, should be developed.

European Union legislation foreshadows the adoption of virus controls for bivalve shellfish when the methods are sufficiently developed [START_REF] Anonymous | European Commission Regulation (EC) No 2073/2005 of 15 November 2005 on microbiological criteria[END_REF] while emergency legislation mandating virus testing on imports to the European Union of strawberries from China [START_REF] Anonymous | European Commission Implementing Regulation (EU)[END_REF] and raspberries from Serbia (Anonymous, 2015) has also been passed in recent years. The availability of a validated EN/ISO standard should enable the introduction of more robust and quality assured food hygiene controls for high risk foods both in Europe and more widely, and will contribute to improved food safety and a reduced burden of food-borne viral illness.

Figure 1

 1 Figure 1 Alignment of GII quantification standard sequences with primer and probes used in the validation study. Matching nucleotides are highlighted in black, mismatches are not highlighted.

Figure 2

 2 Figure 2 Data for part 1; method characterisation for norovirus GI in bottled water. Data points for negative results are shown at 0.01 copies/ml and shaded black. A line of equivalence is shown as a dashed line, and the line of best-fit used to calculate residuals for determination of the LOQ is shown as a solid line.

Figure 3

 3 Figure 3 Probability of detection/limit of detection for norovirus GI in bottled water. The probability of detection function derived from the data is shown as a dotted line. The dashed line marks the 95% probability of detection used to determine the LOD.

Figure 4

 4 Figure 4 Part 2; interlaboratory study results for norovirus GI in bottled water. A, high contamination level; B, medium contamination level; C, low contamination level; D, negative. Quality control results obtained by the expert laboratory are shown to the left of each sub-plot (data points shaded grey). Results for duplicate subsamples as obtained by the 10 participating laboratories are shown to the right (data points unshaded). Data points that were removed from the dataset after quality checking are shaded black. Negative results are shown at 0.01 copies/ml.

Table 1

 1 Expert laboratories involved in the validation

	Laboratory

a nd = not determined

Acknowledgements

The validation of International Standard EN ISO 15216-1 was carried out under the framework of European Mandate No. M381 of DG SANTE and DG GROW (European Commission).

The authors wish to thank the following laboratories for their participation and cooperation in the method characterisation and/or interlaboratory study: Miss N. Gustar, Dr R. Hartnell and Dr D. Lees, Centre for


Guillaume Charria

Sébastien Theetten

Adam Ayouche

Coline Poppeschi

Joël Sudre

Hussein Yahia

Veronique Garcon

François G Schmitt

Sébastien

Multi-scale coastal surface temperature in the Bay of Biscay and the English Channel

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Multi-scale coastal surface temperature in the Bay of Biscay and the English Channel

Guillaume Charria 1 , Sébastien Theetten 1 , Adam Ayouche 1 , Coline Poppeschi 1 , Joël Sudre 2 , Hussein Yahia 3 , Véronique Garçon 2 , François Schmitt The Bay of Biscay and the English Channel, in the North-eastern Atlantic, are considered as a natural laboratory to explore the coastal dynamics at different spatial and temporal scales. In those regions, the coastal circulation is constrained by a complex topography (e.g. varying width of the continental shelf, canyons), river runoffs, strong tides and a seasonally contrasted wind-driven circulation.

Based on different numerical model experiments (from 400m to 4km spatial resolution, from 40 to 100 sigma vertical layers using 3D primitive equation ocean models), different features of the Bay of Biscay and English Channel circulation are assessed and explored. Both spatial (submesoscale and mesoscale) and temporal (from hourly to monthly) scales are considered. Modelled spatial scales, with a specific focus on the variability of fine scale features (e.g. fronts, filaments, eddies), are compared with remotely sensed observations (i.e. Sea Surface Temperature). Different methodologies as singularity and Lyapunov exponents allow describing fine scales features and are applied on both modelled and observed datasets. For temporal scales, in situ high frequency surface temperature measurements from coastal moorings (from COAST-HF observing network) provide a reference for the temporal variability to be modelled. Exploring differences in the temporal scales (from an Empirical Mode Decomposition) advises on the efficiency of our coastal modelling approach.

This result overview in the Bay of Biscay and the English Channel aims illustrating the input of coastal modelling activities in understanding multi-scale interactions (spatial and temporal).