Paolo Ballarini

Centralesupélec France

Analysing oscillatory trends of discrete-state stochastic processes through HASL statistical model checking

The application of formal methods to the analysis of stochastic oscillators has been at the focus of several research works in recent times. In this paper we provide insights on the application of an expressive temporal logic formalism, namely the Hybrid Automata Stochastic Logic (HASL), to that issue. We show how one can take advantage of the expressive power of the HASL logic to define and assess relevant characteristics of (stochastic) oscillators. 1 as opposed to numerical stochastic model checking which requires the complete construction of a model's state-space to assess the exact value of the target measure.

Introduction

Many biological systems are characterised by some form of stochastic oscillatory dynamics. A rough classification may distinguish damped oscillators, that is, systems that exhibit a transient oscillatory trend, from sustained oscillators, that is, systems that exhibit a stationary oscillation lasting in the long run. Assessing the main characteristics of a stochastic oscillator, like, for example, the mean amplitude or the mean period of oscillations, or even the variability of amplitude and period, is not a trivial task.

Oscillation analysis techniques exist both in applied mathematics and signal processing. In mathematics the oscillatory character of systems described in terms of Ordinary Differential Equations Address(es) of author(s) should be given (ODEs) can be assessed through limit-cycle analysis. With signal processing the oscillatory nature of a given signal can be assessed through Fast Fourier Transformation (FFT) or autocorrelation analysis.

In recent times the study of oscillatory systems has attracted the attention of research in the area discrete-state stochastic models, yielding to a number of research works aimed at the application of temporal logic reasoning to characterisation of oscillations [START_REF] Ballarini | Query-based verification of qualitative trends and oscillations in biochemical systems[END_REF][START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF][START_REF] Andrei | Trend-based analysis of a population model of the akap scaffold protein[END_REF][START_REF] David | Runtime verification of biological systems[END_REF]. The goal in that respect is, quite simply, to adapt (stochastic) model checking techniques so that, given a model M , one is capable to obtain answers to questions such as: "does M oscillates?" , "where are the peaks of oscillations ?" "what is the (average) period of oscillations?". Since here we refer to stochastic models answering such questions usually boils down to assessing some distribution of probability (e.g. assessing the steady-state distribution of M , and/or the PDF of the period duration, and/or the PDF of the location of the peaks of oscillations). So far analysis of oscillations through stochastic model checking have been mainly obtained through application of the Continuous Stochastic Logic (CSL) [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] (in some cases joint to its reward-based extensions [START_REF] Kwiatkowska | Stochastic model checking[END_REF]) approach or similarly expressive (linear-time) variants (e.g. Metric Interval Temporal Logic [START_REF] David | Runtime verification of biological systems[END_REF]).

Interestingly, in recent times, Spieler [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF] has shown how CSL model checking can be adapted to oscillation analysis so that both qualitative, as quantitative oscillation related queries can be for-mally assessed against a continuous-time Markov chain (CTMC). This is achieved by coupling of a CTMC model with a timed automaton (TA) "monitor" capable of identifying noisy-periodic traces.

In this paper we extend Spieler's approach by considering a larger class of "monitor" automata, namely that of Hybrid Automata, featured in the Hybrid Automata Stochastic Logic (HASL) [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF], i.e., a recently introduced formalism for assessing a generic class of discrete state stochastic processes, which includes, but is not limited to, CTMCs. Differently from Spieler's approach, which allows one for qualitatively establishing whether a CTMC model is a sustained oscillator and, if so, to assessing some quantitative characteristics (e.g. the average period duration), the methodology we introduce here is limited to assessing quantitative characteristics of a stochastic model which is known (or believed to) oscillates.

Paper contribution. The paper main contribution is one of demonstrating the effectiveness of the HASL formalism as a means to effectively specifying and automatically estimate oscillation related measures. We consider two different approaches: the first is concerned with assessing the oscillation period, the second is concerned with measuring the oscillation peaks (hence the oscillation amplitude). We define two specific types of linear hybrid automata (LHA) that, when synchronised with an oscillatory stochastic process, are capable of detecting the periods, respectively the peaks, of its trajectories and compute on-the-fly classical characteristics like the average duration or the average amplitude of the oscillations as well as more sophisticated ones like the period fluctuation, which allow for assessing the regularity of an oscillator.

We demonstrate such contributions by considering a well-known case study. i.e. the analysis of a model of the circadian clock [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF].

Paper organisation. In Section 2 we introduce the HASL formalism. In Section 3 we describe the basic contribution of the paper, namely the application of HASL to the analysis of oscillations. In Section 4 we demonstrate the HASL-based analysis of oscillations on an example of biological oscillator. We wrap up the paper with some concluding remarks in Section 6.

HASL model checking

HASL framework belongs to the family of so-called statistical model checking methods, whose goal is to produce estimates of a (formally specified) target measure through sampling of the trajectories of the model 1 . HASL is an automata-based type of logic, meaning that it employs automata, and specifically linear hybrid automata (LHA) as machineries for characterising the properties to be investigated. This yields the main feature of HASL, that is, its expressive power which in this paper we are going to demonstrate in respect to the oscillation analysis problem.

Simply speaking the HASL model checking procedure works as follows: given a model D and a certain dynamics of interest encoded by an LHA A, the HASL model checker samples trajectories of the synchronised process D×A, hence selecting only those paths of D that are accepted by A and using them for estimating the confidence-interval of a given target measure (in the following denoted Z), a quantity defined as a function of the LHA variables.

In the remainder of this section we give the basics elements of the HASL formalism, which include: the characterisation of Discrete Event Stochastic Process (DESP), the characterisation of LHA and of the corresponding synchronised process D× A, and that of target measure Z. We refer the reader to [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] for more details.

Discrete Event Stochastic Processes

We refer to a DESP as a discrete-state stochastic process consisting of an enumerable set of states and whose dynamic is triggered by a set of (time consuming) discrete events. We do not consider any restriction on the nature of the distribution associated with events2 . Otherwise said a DESP is a family of random variables {X(t) | t ∈ R ≥0 } where t represents time and where, in the context of this paper, N n is assumed to be the support of X (i.e. we talk in this case of an n-dimensional DESP population model, see Definition 3). Below we give the formal definition of DESP. Such a characterisation is useful to provide us with an algorithmic formulation of the dynamics of a DESP (see below) which is at the basis of the HASL statistical model checking procedure.

Notation. For a generic set A we denote dist(A) the set of possible probability distributions whose support is A, that is, dist(A) = {µ : Σ A → R + |(A, Σ A , µ)} where (A, Σ A , µ) is a probability space for a sigma-algebra Σ A of A. Observe that depending on the nature of A the corresponding probability distributions µ ∈ dist(A) are either continuous (if A is dense) or discrete (if A is finite/discrete).

Definition 1 (DESP) A DESP is a tuple D = S, π 0 , E, Ind, enabled, delay, choice, target where -S is an enumerable (possibly infinite) set of states, -π 0 ∈ dist(S) is the initial distribution on states, -E is a finite set of events, -Ind is a set of functions from S to R called state indicators (including the constant functions), -enabled : S → 2 E are the enabled events in each state with for all s ∈ S, enabled(s) = ∅. -delay : S×E → dist(R + ) is a partial function defined for pairs (s, e) such that s ∈ S and e ∈ enabled(s). -choice : S × 2 E × R + → dist(E) is a partial function defined for tuples (s, E , d) such that E ⊆ enabled(s) and such that the possible outcomes of the corresponding distribution are restricted to e ∈ E . -target : S×E×R + → S is a partial function describing state changes through events defined for tuples (s, e, d) such that e ∈ enabled(s).

A configuration of a DESP consists of a triple (s,τ ,sched) with s being the current state, τ ∈ R + the current time and sched : E → R + ∪{+∞} being the function that describes the occurrence time of each scheduled event (+∞ if an event is not yet scheduled). Observe that a scheduler sched essentially describes the state of the events' queue in a given configuration of the DESP. Thus all (currently) enabled events will have a finite scheduled time sched(e) < ∞, whereas non-enabled events will be associated with an infinite delay, that is, sched(e) = ∞. Thus, within the algorithm for generating a trajectory of a DESP, a scheduler sched provides the occurrence time of the next event to occur (see also Algorithm 1). In the remainder we denote Conf = S × R + × Sched the set of possible configurations of a DESP (where Sched denotes the set of possible schedules functions for the events of the DESP). Also for a configuration c = (s , τ , sched) ∈ Conf , we denote c(s), c(τ ) and c(sched) the state s, respectively the time τ and the schedule sched of configuration c.

For a state s, enabled(s) is the set of events enabled in s. For e ∈ enabled(s), delay(s, e) is the distribution of the delay between the enabling of e and its possible occurrence. Furthermore, if we denote δ m the delay of the earliest event in the current configuration (s,τ ,sched) of the process, and E min ⊆ enabled(s) the set of events with earliest delay, then choice(s, E min , δ m ) describes how the conflict between the concurrent events in E min is randomly resolved: i.e. choice(s, E min , δ m )(e ) is the probability that e ∈ E min will be selected hence occurring with delay δ m . Finally function target(s, e, d) denotes the target state reached from s on occurrence of e after waiting for d time units.

Dynamics of a DESP. The evolution of a DESP D can be informally summarised by an iterative procedure consisting of the following steps (assuming (s , τ , sched) is the current configuration of D): 1) determine the set E min of events enabled in state s and with minimal delay δ m ; 2) select the next event to occur e next ∈ E min by resolving conflicts (if any) between concurrent events through probabilistic choice according to choice(s, E min , τ ); 3) determine the new configuration of the process re-sulting from the occurrence of e next , this in turns consists of three sub-steps: 3a) determine the new state resulting from occurrence of e next , i.e. s = target(s, e next , δ m ); 3b) update the current time to account for the delay of occurrence of e next , i.e. τ = τ + δ m ; 3c) update the schedule of events according to the newly entered state s (this implies setting the schedule of no longer enabled events to +∞ as well as determining the schedule of newly enabled events by sampling through the corresponding distribution). Such procedure is (semi-formally) summarised in Algorithm 1. A path (or trajectory) of a DESP is a sequence of configurations σ = c 1 , c 2 , c 3 , . . . resulting from the execution of the procedure highlighted by Algorithm 1. We formalise this in the following definition. The notion of DESP path will be used later on for reasoning about the dynamics of a DESP and in particular for reasoning about oscillations.

Definition 2 (Path of a DESP) For a DESP D = S, π 0 , E, Ind, enabled, delay, choice, target we define the set of finite paths as P ath * ⊆ n∈N Conf n , with Conf the set of configurations of D. We denote σ = (c 0 , c 1 . . . , c n ) ∈ P ath * , where π 0 (c 0 (s)) > 0 and

∀0 ≤ i < n, ∃e ∈ E min (c i (s)) such that c i+1 (s) = target(c i (s), e, c i (τ ))
. By extension we denote P ath ω as the set of infinite path and P ath = P ath * ∪ P ath ω as the of all paths of a DESP.

In the remainder we might refer to a DESP path using σ = (c 0 , c 1 . . . , c n ) or, depending on the context, simply indicating the corresponding sequence of states σ = (c 0 (s), c 1 (s) . . . , c n (s)), or simply σ = (s 0 , s 1 . . . , s n ). Furthermore for a path σ = (c 0 , c 1 . . . , c n ) we use the following notations: for i ∈ N, σ[i] = c i (s) denotes the i-th state, time(σ, i) denotes the time spent in state σ[i], while for t ∈ R + , σ@t denotes the state in which σ is at time t, that is, σ@t = σ[j] such that j is the smallest j with t ≤ i≤j i=0 time(σ, i). Since in this paper we deal with the analysis of discrete-state biological models representing the evolution of the molecular population of n species, we introduce the notion of DESP population model.

= (s 1 , . . . , s i , . . . s n ) ∈ S of D corresponds to state s i ∈ S i of D i .
By extension for σ ∈ P ath a path of a DESP population model we denote σ i the i th projection of σ, thus if σ = (s 1 0 , . . . s n 0 ), (s DESP in terms of GSPN For implementation convenience, in the context of HASL and in particular of the associated model checking tool COSMOS [START_REF]COSMOS home[END_REF][START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF], we represent DESP models in terms of stochastic petri nets, and more precisely we adopt (the non-markovian extension3 of) Generalised Stochastic Petri Net (GSPN) [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] as the high-level input formalism for expressing a DESP model. Thus, in this context, DESP indicators are actually GSPN indicators, that is: they are expressions which contain references to the (marking of the) places of a GSPN model. For the sake of brevity here we assume familiarity with the GSPN formalism, referring the reader to the literature [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] for details. GSPN semantics is briefly presented later on through description of a simple GSPN model (see Figure 1).

Example: DESP indicators within LHA. In the LHA of Figure 1 (right) the indicator protA, which refers to the marking of the GSPN place named protA in Figure 1 (left), is used within the updates of the self-loop edges of location l 0 . Specifically indicator protA is used to update the LHA variable a with the current number of tokens contained in GSPN place protA. In the following we formally introduce the notion of synchronised LHA and then informally describe the stochastic process resulting from the product of a DESP and a synchronised LHA.

Definition 5 A synchronised linear hybrid automaton is a tuple A = E, L, Λ, I, F, X, flow, → where:

-E is a finite alphabet of events; -L is a finite set of locations; -Λ : L → P rop is a location labelling function; -I ⊆ L is the initial locations; -F ⊆ L is the final locations; -X = (x 1 , ...x n ) is a n-tuple of data variables; flow : L → Ind n associates an n-tuple of indicators with each location (the i th projection flow i denotes the flow of change of variable

x i ). -→⊆ L × (2 E × Const) ({ } × lConst) × Up × L is the set of edges of the LHA ,
where denotes the disjoint union, Const and lConst denotes the set of possible constraints, respectively left closed constraints, associated with A (see details below), Up is the set of possible updates for the variables of A and P rop ⊆ Ind denotes the subset of boolean valued DESP indicators.

Before presenting informally the synchronisation of a DESP with an LHA we start by describing the various parts of an LHA. In what follows we denote indicators symbolically by greek letters α, α ∈ Ind, while we use capital letters A, B, . . . to refer to names of GSPN places (within concrete indicators instances) and x 1 , x 2 , . . . to denote LHA variables. Thus, for example, α ≡ A + 2B is an indicator whose value is given by the sum of the marking of place A with twice the marking of place B.

Location proposition: function Λ associates each location l ∈ L with a proposition (also called location invariant in the remainder) Λ(l) ∈ P rop representing a condition under which a location can be entered. A location proposition consists of a boolean combination of inequalities involving DESP indicators and has the following form Λ(l) ≡ i (α i ≺ α i ) with α i , α i ∈ Ind, and ≺∈ {=, <, >, ≤, ≥}. Notice that indicators can be constant functions, thus, for example, a location proposition may consist of comparing indicators' values against constant thresholds, as in, e.g., Λ(l) ≡ A ≥ 10, or it may consist of comparing different indicators one another, as in, e.g., Λ(l)

≡ A ≤ B or Λ(l) ≡ A ≤ B 2 √ C.
A location proposition (given by Λ) is shown by a label next to the location it refers to. For convenience no label is shown next to unconstrained locations, i.e., locations associated to a tautology like ≡ (α i = α i ). Location propositions are evaluated against states of a DESP. Thus for s ∈ S a DESP state and l ∈ L a location of an LHA we say that s satisfies the invariant Λ(l), denoted s |= Λ(l), if Λ(l)(s) = true (where Λ(l) is the value of the boolean expression obtained by replacing each indicator α ∈ Λ(l) with its value α(s)). Furthermore given two edge locations l and l we say that the their respective invariants are inconsistent, denoted Λ(l) ∧ Λ(l ) ⇔ false, if there cannot exist a state s that satisfies Λ(l) ∧ Λ(l ). For example, if Λ(l) ≡ A ≤ 2 and Λ(l ) ≡ A > 2 then trivially A ≤ 2 ∧ A > 2 ⇔ false. This means that l and l are mutually exclusive, which is a necessary condition for LHA with multiple initial locations (see conditions c1 below).

Edge constraint: Edge constraints describe necessary conditions for an edge to be traversed. We denote Const (resp. lConst) the set of constraints (resp. left-closed constraints) of an LHA edge. An edge constraint consists of a boolean combination of inequalities involving both DESP indicators and LHA variables. They have the following form γ ≡ j ( 1≤i≤n α ij x i ≺ α j ) with α ij , α j ∈ Ind, x i ∈ X and ≺∈ {=, <, >, ≤, ≥}. Simple examples of edge constraints can be: γ ≡ (2x 1 + 3x 2 ≤ 5) or also γ ≡ (Ax 1 = 5), where A in the latter case represents a DESP indicator (i.e., the marking of the GSPN place named A). Given a location l of the LHA and a state s of the DESP, the inequalities γ j ≡ 1≤i≤n α ij x i ≺ α j of constraint γ, evolve linearly with time. Thus each γ j gives an interval of time during which γ j is satisfied, and the complete constraint γ is satisfied over the intersection of the intervals satisfying the conjuncts γ j . We say that a constraint is left closed if, whatever the current state s (defining the values of the DESP indicators), the time at which the constraint is satisfied is a union of left closed intervals (for example, γ ≡ (x 1 ≥ 5) is left-closed whereas γ ≡ (x 1 > 5) is not). We denote lConst ⊆ Const the subset of left-closed constraint. For efficiency the constraint of autonomous-edges (see below) must be left-closed. Edge constraints are evaluated against pairs (s, ν) ∈ S ×V al where s ∈ S is a state of a DESP and ν : X → R ∈ V al is a valuation that maps every LHA data variable to a real value (we denote V al the set of all possible valuations). For ν ∈ V al, ν(x) denotes the value of variable x through valuation ν. Given γ j ≡ 1≤i≤n α ij x i ≺ α j an inequality contained in an edge-constraint γ ≡ j γ j , its interpretation w.r.t. ν and s, denoted γ j (s, ν), is defined by γ j (s, ν) = 1≤i≤n α ij (s)ν(x i ) ≺ α j (s). We write (s, ν) |= γ j if γ j (s, ν) = true and, by extension, (s, ν) |= γ iff (s, ν) |= γ j for all j. Furthermore given two edge constraints γ and γ we say that their conjunction γ ∧γ is inconsistent, denoted γ∧γ ⇔ false, if there exists no combination (s, ν) ∈ S ×V al that satisfies it. For example, if γ ≡ x 1 ≤ 2 and γ ≡ x 2 > 2 are the constraints for two edges then trivially (x 1 ≤ 2) ∧ (x 1 > 2) ⇔ false, meaning the two edges cannot be concurrently enabled (see conditions c2 and c3 below).

Edges update: an edge update U = (u 1 , ..., u n ) ∈ U p is an n-tuple of functions characterising how each LHA variable x k is going to be updated on traversal of the edge. Each function u k (1 ≤ k ≤ n) of an edge update U = (u 1 , ..., u n ) ∈ U p is of the form x k = 1≤i≤n α i x i + c where the α i and c are DESP indicators. Similarly to edge constraints, updates are evaluated against pairs (s, ν) ∈ S×V al. Given an update U = (u 1 , . . . , u n ), we denote by U (s, ν) the valuation defined by U (s, ν)

(x k ) = u k (s, ν) for 1 ≤ k ≤ n.
Locations flow: a location flow is an n-tuple of indicators f low(l) = (α 1 , . . . , α n ), where α i ∈ Ind describes the gradient at which variable x i ∈ X changes while the automaton sojourns in location l. Specifically when location l is entered the rate of change of each x i is established by the valuation, w.r.t. to the state the DESP is at on entering of l, of the corresponding α i . Observe that, if each α i in f low(l) is a constant function (e.g. α i = c i , with c i ∈ R) then each variable x i changes at constant rate throughout the sojourn in l. However this is not necessarily the case for variables whose flow is given by a non-constant indicator, like, for example, α i = c 1 A + c 2 B, with c 1 , c 2 ∈ R and A, B representing the marking of two GSPN places named A and B. In this case the flow of change of x i depends on the marking of places A and B, and such marking may change during the sojourn in l, for example if a synchronising self-loop edge l → l exists which synchronises with some DESP event whose occurrence modify the marking of A or B.

Having described the DESP indicators dependent elements of an LHA we now see how they are all combined within the characterisation of an LHA edge.

Edges of an LHA. An edge l E ,γ,U ----→ l of an LHA is labelled by: 1) a constraint γ, 2) a set of event labels E , 3) an update U . An edge can be either synchronous or autonomous. A synchronous edge is one whose traversal is triggered by the occurrence of an event of the DESP in particular an event e ∈ E ⊆ E where E is the set of event names labeling the edge. An autonomous edge, on the other hand, is one whose traversal is independent of the occurrence of DESP events, hence the event label for autonomous edges is E ≡ , where is the label used for representing a "pseudo-event".

The class of LHA for HASL is further restrained by the following conditions:

-c1 (initial determinism): ∀l = l ∈ I, Λ(l) ∧ Λ(l ) ⇔ false. This means that independently of the interpretation of the indicators, hence of the synchronising DESP model, at most one initial location l ∈ I can have its constraint

Λ(l) verified. -c2 (determinism on events:) ∀E 1 , E 2 ⊆ E : E 1 ∩ E 2 = ∅, ∀l, l , l ∈ L, if l E1,γ,U ----→ l and l E2,γ ,U -----→ l are two distinct transi- tions, then either Λ(l) ∧ Λ(l ) ⇔ false or γ ∧ γ ⇔ false.
Again this equivalence must hold whatever the interpretation of the indicators occurring in Λ(l), Λ(l ), γ and γ .

-c3 (Determinism on :) ∀l, l , l ∈ L, if l ,γ,U ---→ l and l ,γ ,U ----→ l are two distinct transitions, then either Λ(l) ∧ Λ(l ) ⇔ false or γ ∧ γ ⇔ false.
-c4 (no -labelled loops:) For all sequences

l 0 E0,γ0,U0 ------→ l 1 E1,γ1,U1 ------→ • • • En-1,γn-1,Un-1 -----------→ l n such that l 0 = l n , there exists i ≤ n such that E i = .
Synchronisation of LHA and DESP. The role of a synchronised LHA A is to select specific trajectories of a corresponding DESP D while collecting relevant data (maintained in the LHA variables) along the execution. Synchronisation is technically achieved through the product process D×A whose formal characterisation, for the sake of brevity, we omit in this paper: we provide however an intuitive description of the D × A semantics.

The product D × A is itself a DESP whose states are triples (s, l, ν) where s is the current state of the D, l the current location of the A and ν : X → R the current valuation of the variables of A. Formally the set of states of the product process D ×A is defined as S = (S × L × V al) {⊥}, where V al denotes the set of possible variables' valuations and ⊥ denotes the rejecting state, i.e., the state entered when synchronisation fails, hence when a trajectory is rejected (see below). Notice that a configuration of the product DESP D × A has the following form ((s, l, ν), τ, sched ), where (s, l, ν) is the current state of D × A, τ ∈ R + is the current time, and sched is the schedule of the enabled events of D × A. The synchronisation starts from the initial state (s, l, ν), where s an the initial state of the DESP (i.e. π 0 (s) > 0), l is an initial location of the LHA (i.e. l ∈ I) and the LHA variables are all initial set to zero (i.e. ν = 0) 4 .

From the initial state the synchronisation process evolves through transitions where each transition corresponds to traversal of either a synchronised or an autonomous edge of the LHA 5 . Furthermore if an autonomous and a synchronised edge are concurrently enabled the autonomous transition is taken first. Let us suppose that (s, l, ν) is the current state of process D×A and describe how the synchronisation evolves. If in the current location of the LHA (i.e. location l of the current state (s, l, ν)) there exists an enabled autonomousedge l ,γ,U ---→ l , then that edge will be traversed leading to a new state (s, l , ν ) where the DESP state (s) is unchanged whereas the new location l and the new variables' valuation ν might differ from l, respectively ν, as a consequence of the edge traversal. On the other hand if an event e of process D (corresponding to transition s e -→ s of D) occurs in state (s, l, ν), either an enabled synchronous edge l E ,γ,U ----→ l (with e ∈ E ) exists leading to new state (s , l , ν ) of process D×A (from which synchronisation will continue) or the synchronisation halts hence the trace is rejected (formally this is achieved with the system entering the rejecting state ⊥).

Enabling of an LHA edge. Let us briefly describe how the enabling, hence the traversal, of an LHA edge is established. Let (s, l, ν) be the current state of process D × A. An edge l E,γ,U ----→ l being it autonomous or synchronous originating in l is en-abled if the following two conditions hold: 1) if the edge constraint is satisfied in state (s, l, ν) (i.e., if (s, ν) |= γ) 2) if the location invariant of the target location Λ(l ) is satisfied in the state s reached by traversal of the edge, i.e., if s |= Λ(l ) (observe that if the considered edge is autonomous then necessarily s = s, whereas if it is synchronous then possibly s = s). Finally for a synchronous edge to be enabled, in addition to 1) and 2), it must be the case that the DESP event e occurring while in (s, l, ν) is captured by the the edge, i.e., e ∈ E.

Remarks. The above described synchronisation of a DESP and an LHA, which HASL model checking is based on, requires certain properties to hold, namely: uniqueness, convergence and termination of the synchronisation. This means that for A a synchronised LHA then for any (infinite) path σ of a synchronising DESP model: 1) there must be exactly one synchronisation with A, 2) synchronisation cannot go on indefinitely due to an infinity of consecutive autonomous events, 3) path σ should lead to an absorbing state (i.e. a final location of the A or the rejecting state ⊥) with probability 1. The uniqueness property is guaranteed by constraint c1, c2 and c3 of the LHA definition whereas convergence is a consequence of constraint c4. On the other hand termination of the synchronisation is not explicitly guaranteed, however can be ensured by structural properties of A and/or D.

HASL expressions. The second component of an HASL formula is an expression related to the automaton. Such an expression, denoted Z, is defined by a specific grammar [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] of which here we consider only the basic elements given in [START_REF]COSMOS home[END_REF].

Z ::= E[Y ] | P Y ::= last(y) | min(y) | max(y) | avg(y) y ::= c | x | y + y | y × y | y/y (1)

Z is either either an expectation expression E[Y ], or a probability expression P . An expectation expression Z = E[Y ] represents the expected value

Example (synchronisation of DESP and LHA): To understand how synchronisation of a DESP with an LHA works let us consider a simple example. Figure 1 depicts a toy DESP model in GSPN form (on the left) coupled with a simple LHA (on the right). The GSPN model represents the basic steps of gene expression: 1) binding/unbinding of an activator protein to the promoter of gene A; 2) of transcription of a gene into an mRNA molecule; 3) degradation of the mRNA 4) translation of the mRNA into the expressed protein prot A. The states of the DESP consist of 4-tuples s = (protA, geneA ,A geneA, mrnA) ∈ N 4 , corresponding to the marking of the 4 places of the GSPN, whereas the event set is E = {bind, unbind, degrade, transc, transl}, corresponding to the 5 timed-transitions of the GSPN. The LHA A, on the other hand, consists of: two locations, l 0 (initial) and l 1 (final), and three data variables t (a clock), n (for counting the number of occurrences of the transc event) and a (for keeping track of the population of prot A). Notice that the invariant of both locations is Λ(l 0 ) = Λ(l 1 ) = , (hence no label is associated to l 0 , l 1 ), meaning that both locations can be entered without constraint. The initial state of the product process D × A is s 0 = ((2, 1, 0, 0), l 0 , ν 0 ), where ν 0 is the zero

valuation (i.e., ν 0 (t) = ν 0 (n) = ν 0 (a) = 0). A has two synchronised (self-loop) edges l 0 {transc},n<N,{n++,a=protA} ---------------------→ l 0 ,
which synchronises with occurrences of the transc event, and l 0 E\{transc},n<N,{a=protA} -------------------→ l 0 , which synchronises with occurrences of any other event but transc (i.e., E \ {transc}). The constraint for both synchronised edges is n < N which means they can be traversed as long as the number of observed occurrences of transc, which is stored in n, is less than N . Both updates for the two synchronised edges refer to a single indicator, namely protA, whose value is given by the marking of the GSPN place labelled protA, but they are slightly different. The update for the edge which synchronises with transc is {n ++, a = protA}, meaning that whenever the edge is traversed (i.e., on occurrence of a transc event) the counter n is incremented and the current marking of place protA is stored in a. On the other hand the update for the edge which synchronises with E \ {transc} the update is simply {a = protA} as clearly n must be incremented only on occurrence of transc. Furthermore A has an autonomous edge l 0 ,(n=N ),∅ ------→ l 1 leading to the final location l 1 . Such edge gets enabled as soon as its constraint (n = N ) is satisfied, that is, as soon as a state s N = ((n 1 , n 2 , n 3 , n 4 ), l 0 , ν N ), is reached with ν N being a valuation such that ν N (n) = N . In any such state s N the autonomous edge is traversed (leading to state sstop = ((n 1 , n 2 , n 3 , n 4 ), l 1 , ν N )) and the synchronisation stops. of a random variable Y built on top of basic path operators (last(y), min(y), max(y) , avg(y)). Each such path operator take as argument y an algebraic combination of the LHA data variables x, and is evaluated along a (synchronised) path that is accepted by the automaton. Intuitively the meaning of path operators is as follows: last(y) represents the value that expression y has at the instant a path is accepted, while min(y) (max(y), respectively avg(y)) represents the minimum (maximum, respectively average) value assumed by y along an accepted path. Expression Z = P , on the other hand, simply represents the probability that a path is accepted by the LHA. This is given by the ratio between the number of accepted paths and total number of paths generated throughout a simulation experiment.

In recent updates the COSMOS model checker [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF] has been enriched with facilities for assessing the Probability (Cumulative) Distribution Function (PDF, respectively CDF) of the value that an expression Y takes at the end of a synchronising path. Notice that PDF and CDF HASL expressions, are only high-level macros supported by the COSMOS tool in order to give the user the possibility to specify PDF/CDF measures6 in a straightforward manner. Thus COSMOS supports the following syntax for estimating a PDF measure: Z = P DF (Y, s, l, h), where Y is the path dependent expression whose PDF is to be estimated while l and h are the lower, respectively higher, bound of the interval representing the support of Y (i.e. estimation of the PDF of Y is done assuming Y takes value in [l, y]) and s < (hl) is the width of each sub-interval in which the considered support [l, y] is discretised. Thus during estimation of Z = P DF (Y, s, l, h) COSMOS internally maintains a counter for each of the (hl)/s sub-intervals. Each such counter is incremented if the value of Y on acceptance of a trace falls in the corresponding sub-interval. Then the value returned by COSMOS for Z = P DF (Y, s, l, h) is the array of frequencies obtained by dividing each of the above counters by the total number of generated trajectories.

Example. We consider some examples of HASL formula referred to the GSPN model of Figure 1 (left). Observe that all three formulae φ 1 , φ 2 , φ 3 use the same monitor automaton, i.e., A (Figure 1, right). This means that φ 1 , φ 2 , φ 3 will all select the same type of trajectories, specifically those trajectories whose prefix contains exactly N occurrences of the transc event.

For φ 1 the expression to be estimated is

Z 1 = E[last(t)],
which is, the expected value that the LHA variable t exhibits at the end (last(t)) of an accepted, synchronised trace. Thus for each trace σ sampled from process D × A the value of t at the moment σ is accepted (i.e., on occurrence of the N -th transc event) is retained as a sample for the confidence-interval estimation of the expectation of last(t). For φ 2 the expression to be estimated is Z 2 = E[max(a)], which is, the expected value of the maximum that LHA variable a exhibited along a synchronised trace. The maximum of an LHA variable along a trace is automatically computed on-the-fly during the sampling of a trace hence the value max(a) for a sampled trace σ is known straight away on acceptance of σ. Notice that a stores the marking of the GSPN place prot A, therefore Z 2 represents the expected value of the maximum number of protein A observed over the first N transcription events. Finally for φ 3 the expression Z 3 = P DF (last(t), 0.1, 0, 10) corresponds to estimating with what probability the value of last(t) (i.e., the value of t at the end of a sampled trajectory) falls within a discretised subinterval of [0, 10]. In this case we consider k = (10 -0)/0.1 = 100 sub-intervals of ∆ = [0, 10] each of width 0.1 and with the k-th subinterval being

∆ k = [0+0.1•k, 0+0.1•(k+1)] with 0 ≤ k ≤ 99.
In practice, for estimating PDF-expressions like Z 3 , COSMOS uses k internal variables, which we may call N ∆ k last(t) each of which counts how many times the value of last(t) observed at the end of a sampled trajectory σ has been found falling into the k-th interval ∆ k . The probability that last(t) ∈ ∆ k then simply corresponds to N ∆ k last(t) /n, where n is the total number of sampled trajectories. Thus the estimate of Z 3 produced by COSMOS is the k-tuple of variables (N ∆0 last(t) /n, . . . , N ∆ k last(t) /n, . . . N ∆99 last(t) /n).

COSMOS statistical model checker

COSMOS 7 [6] is a prototype software platform for HASL-based statistical model checking. It employs confidence interval techniques for estimating the mean value of relevant performance measures expressed in terms of HASL formulae against a given GSPN model. COSMOS has been recently integrated in the CosyVerif platform [START_REF]CosyVerif home[END_REF] which adds to the original command line interface (available with the first version) the possibility of drawing the input elements (i.e. GSPN and LHA) through a user a graphical interface. Software platforms featuring statistical model checking functionalities similar to COSMOS include: PRISM [START_REF]Prism home page[END_REF], UPPAAL-SMC [START_REF] Bulychev | Uppaalsmc: Statistical model checking for priced timed automata[END_REF], and PLASMA [START_REF] Jégourel | A platform for high performance statistical model checking -plasma[END_REF], APMC [START_REF] Herault | APMC 3.0: Approximate verification of discrete and continuous time Markov chains[END_REF], YMER [START_REF] Younes | Ymer: A statistical model checker[END_REF], MRMC [START_REF] Katoen | Simulation-based CTMC model checking: an empirical evaluation[END_REF] and VESTA [START_REF] Sen | VESTA: A statistical model-checker and analyzer for probabilistic systems[END_REF]. We refer the reader to [START_REF]COSMOS home[END_REF][START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF] for more details on COSMOS.

Measuring oscillations with HASL

Intuitively an oscillation is the periodic variation of a quantity around a given value. In mathematical terms this is associated with the definition of (non-constant) periodic function. i.e. function f :

R + → R for which ∃P ∈ R + such that ∀t ∈ R + , f (t) = f (t+P )
, where P is called the period (e.g. trace in Figure 2(a)). In the context of stochastic models such a "deterministic" characterisation of periodicity is of little relevance, as the trajectories of a stochastic oscillator being strictly periodic (as in f (t) = f (t + P )), will have (unless in degenerative cases) zero probability. More generally the traces of (discrete-state) stochastic oscillators are characterised by a remarkable level of noise (e.g. trace in Figure 2(b)). Generally speaking, oscillation can either be a transient behaviour (a model which oscillates for a finite duration) or a limiting behaviour (i.e. a model that exhibits a sustained oscillation). In [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF] Spieler, whose work tackles CSL based analysis of CTMC oscillators, characterised sustained oscillations as the absence of both divergence and convergence, meaning that a stochastic model that oscillates sustainably is one whose trajectories σ

do not diverge (lim t→∞ σ(t) < ∞) nor converge ( n ∈ N : lim t→∞ σ(t) = n).
Otherwise said a stochastic model oscillates sustainably if and only if the probability measure of the converging trajectories and diverging trajectories is null [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF].

In order to study the dynamics of stochastic oscillators, in the following we introduce two (orthogonal) characterisations of oscillatory traces. The first one, which following [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF] we call noisy periodicity, allows us for observing the period duration of an oscillator. The second allows us to locate the maximal and minimal peaks of oscillating traces. We first recall the definition of trajectory of a DESP.

Measuring the period of oscillations

As we previously argued, the mathematical characterisation of periodic function is of little use when it comes with analysing the periodic nature of a stochastic model. As a consequence here we consider an alternative characterisation of periodicity which is suitable for capturing the noisy nature of stochastic oscillations. For this we establish a partition of a DESP state-space induced by two threshold levels L, H ∈ N with L < H and we say that, with respect to a specific observed species (i.e. one of the n dimensions of the DESP) a trajectory oscillates or, equivalently is noisy periodic, if it traverse Definition 6 (noisy periodic trajectory) For D an n-dimensional DESP population model, let dom(D) be its state space, and dom i (D) the projection of dom(D) along the i th dimension 1 ≤ i ≤ n of D. Let L, H ∈ N, L < H, be two levels establishing the partition dom i (D) = low ∪ mid ∪ high with low = (-∞, L), mid = [L, H) and high = [H, ∞). A trajectory σ of D is said noisy periodic with respect to the i th dimension, and with respect to the considered L, H induced partition of dom i (D) if the projection σ i visits the intervals low), mid and high infinitely often.

In the remainder rather than referring to the periodicity with respect to the i th dimension we refer to the periodicity with respect to the population of species A, where A is the symbolic name of the observed species corresponding to one of the Petrinet place in the GSPN representation of D. Thus, with a slight abuse of notation, we will denote σ A a trace which is noisy periodic w.r.t. species A.

Given a noisy periodic trace we are interested in measuring the basic characteristics of its oscillatory nature, such as, the (average) duration of the oscillation period. For this we first need to establish what we mean by period. Intuitively a period, for a trace which is noisy periodic (in the sense of Definition 6), corresponds to the time interval between two consecutive sojourns in one of the two extreme regions of the partition (e.g., low region), interleaved by (at least) one sojourn into the opposite region (e.g., high region). Figure 3 gives an example of period realisations for a noisy periodic trace: the first two period realisations, denoted p1 and p2, are delimited by the mid-to-low crossing points corresponding to the first entering of the low region which follows a previous sojourn in the high region. Such an intuitive description of noisy period of a noisy periodic trace is formalised in Definition 7. We first introduce the notion of crossing points sets associated to a noisy periodic trace.

Given a noisy periodic trace σ A we denote τ j↓ (respectively τ j↑ ), the instant of time when σ A enters for the j-th time the low (respectively the high) region. T ↓ = ∪ j τ j↓ (resp. T ↑ = ∪ j τ j↑ ) is the set of all low-crossing points (reps. high-crossing points). Observe that T ↓ and T ↑ reciprocally induce a partition on each other. Specifically T ↓ = ∪ k T k↓ where T k↓ is the subset of T ↓ containing the k-th sequence of contiguous low-crossing points not interleaved by any high-crossing point. Formally

T k↓ = {τ i↓ , . . . , τ (i+h)↓ |∃k , τ (i-1)↓ < τ k ↑ < τ i↓ , τ (i+h)↓ < τ (k +1)↑ }. Similarly T ↑ is partitioned T ↑ = ∪ k T k↑
where T k↑ is the subset of T ↑ containing the k-th sequence of contiguous high-crossing points not interleaved by any low-crossing point. For example, with respect to trace σ A depicted in Figure 3 we have that trace can be seen as a collection of realisations of certain random variables. For example the instants of time τ j↓ , τ j↑ are realisations of the random variables (which we could denote x j↓ , respectively x j↑ ) corresponding to the timing of entering the low, respectively high, regions. Similarly the duration of the k-th period contained in a trace can be seen as the realisation of a random variable 8 . We formalise the notion of noisy period realisation in the next definition.

T ↓ = T 1↓ ∪ T 2↓ ∪ T 3↓ . . . with T 1↓ = {τ 1↓ , τ 2↓ }, T 2↓ = {τ 3↓ }, T 3↓ = {τ 4↓ }, while T ↑ = T 1↑ ∪ T 2↑ ∪ T 3↑ . . . with T 1↑ = {τ 1↑ }, T 2↑ = {τ 2↑ }, T 3↑ = {τ 3↑ }. Observe that a noisy periodic p1 L H low mid high p2 p0 T ⌧ 2" ⌧ 3" ⌧ 1" ⌧ 1# ⌧ 2# ⌧ 3# ⌧ 4# A Fig. 3
Definition 7 (k th noisy period realisation) For σ A a noisy periodic trajectory with crossing point times T ↓ = ∪ k≥1 T k↓ , respectively T ↑ = ∪ k≥1 T k↑ , the realisation of the k th noisy period, denoted t p k , is defined as t p k = min(T (k+1)↓ )min(T k↓ ).

Observe that in Definition 7 t p k could alternatively be defined as t p k = min(T (k+1)↑ )-min(T k↑ ), that is, with respect to crossing into the high region, rather that into the low region. It is straightforward to show that both definitions are semantically equivalent, i.e., the average value of t p k measured along a trace (see Definition 8) is equivalent with both definitions. Furthermore observe that a noisy periodic trace (as of Definition 6) contains infinitely many realisations of (noisy) periods. In the remainder we will refer to the N -prefix of a noisy periodic trace σ A , meaning the prefix of σ A that consists of the first N noisy period realisations.

As an example of period realisations, let us consider the noisy periodic trace in Figure 3 whose first two period realisations are t p1 = τ 3↓τ 1↓ and t p2 = τ 4↓τ 3↓ . Notice that the time interval denoted as p0 in Figure 3 does not represent a complete period realisation as there's no guarantee that T = 0 corresponds with the actual entering into the low region. Definition 7 correctly does not account for the first spurious period p0.

Having introduced the notion of noisy period realisation we now look at the problem of estimating two characteristic measures related to it, namely, the period average and the period fluctuation. By period average we simply mean the average value of the period realisations sampled along a trace. On the other hand by period fluctuation we mean a measure of the variability of the period realisations along a trace. Observe that, from the point of view of analysis, period fluctuation allows us to analyse the regularity of the observed oscillator. In this respect a "regular" oscillator is one whose traces consists of little variable periods (i.e., small fluctuation), as opposed to an "irregular" one whose traces exhibits variable periods (i.e., large fluctuation). We demonstrate the analysis of oscillation regularity through fluctuation assessment in Section 4.

Definition 8 (period average) For σ A a noisy periodic trajectory the period average of the first n ∈ N period realisations, denoted t p (n), is defined as

t p (n) = 1 n n k=1 t p k , where t p k is the k-th period realisation.
Observe that for a sustained oscillator, the average value of the noisy-period, in the long run, corresponds to the limit t p = lim n→∞ t p (n).

Definition 9 (period fluctuation) For σ A a noisy periodic trajectory the period fluctuation of the first n ∈ N period realisations, denoted

s 2 tp (n), is de- fined as s 2 tp (n) = 1 n n k=1 (t p k -t p (n)) 2
, where t p k is the k-th period realisation and t p (n) is the period average for the first n period realisations.

Note that the period fluctuation is in essence defined as the variance of the period realisations along a trace. In the remainder we show how, through automaton A per , we can estimate the period fluctuation on-the-fly, that is, as the noisy periodic traces are generated and scanned by A per . For this we employ an adaptation of the so-called online algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF] for computing the variance out of a sample of observations.

In the following we introduce an LHA automaton, called A per , which is targeted to estimating both the average and the fluctuation of the first N the period realisations occurring along the simulated noisy periodic traces.

The automaton A per . The LHA A per depicted in Figure 4 is designed for detecting noisy periods of an observed species (here denoted A). It consists of an initial transient filter (locations l 0 , l 0 ) plus three main locations low, mid and high (corresponding to the partition of A's domain induced by thresholds L < H). The intuition behind the structure of the A per automaton is as follows: the transient filter is used to simply let the simulated ----→ l 0 , simply observing the occurrences of any event of E (i.e., the event set of synchronised GSPN-DESP model), hence letting the simulated trajectory unfold for a fixed duration given by parameter initT : when t ≥ initT A per moves, through autonomous edges, to either l 0 , if by t = initT the simulated trajectory is not in a state of the low = (-∞, L] region (i.e., if the invariant A > L of l 0 is satisfied), or to location low if the current state of the trajectory belongs to the low = (-∞, L] region (i.e., if the invariant A ≤ L of low is fulfilled). If l 0 is entered then the simulated trace is let further unfolding (synchronised self-loop l 0 E, ,∅ ----→ l 0 ) until a state within low = (-∞, L] is reached, in which case the invariant of location low is fulfilled hence the autonomous edge l 0 , ,...

----→low is traversed. Observe that on entering of low the global timer variable t is reset and the period counter n is initialised to -1 (this is so to avoid the first spurious period, denoted p0 in Figure 3, to be considered amongst the detected ones). Once in location low the actual detection of the period realisations begins 9 and the automaton gets looping between the low, mid and high locations for as long as N periods have been detected. From low the automaton follows the profile exhibited by the ob- 1 The data variables of automata Aper of Figure 4 for measures of noisy-periodicity served population A, thus moving to mid (and possibly back) as soon as the population of A grows and a state of the mid = (L, H) region is entered (i.e., corresponding to the L < A < H invariant of mid location becoming satisfied), and then to high (and possibly back) as soon as the population of A enters the high = [H, +∞) region (corresponding to the A ≥ H invariant of high location). On entering the high location the boolean variable top is set to true (i.e., top = 1). This allows then for distinguishing between the mid-to-low transitions of kind mid E,(...∧top=1),...

----------→low, which correspond to an actual closure of a period realisation (i.e., those τ j↓ preceded by a sojourn in the high = [H, +∞) region), from those of kind mid E,(...∧top=0),...

----------→low

which correspond to a return to low without having previously sojourned in high. Observe that from mid location there are four possible (mutually exclusive) ways of entering the low location. If the sojourn in mid has not been preceded by a sojourn in high edge mid E,(n<N ∧top=0),...

------------→low is enabled. On the other hand if the sojourn in mid has been preceded by a sojourn in high but low is going to be re-entered for the first time (i.e., n = -1) then the timer t is reset (representing the start time of actual period detection) and the counter of detected periods n is set to zero (again representing the actual beginning of counting of period detection). Furthermore if the sojourn in mid has been preceded by a sojourn in high and the period to be detected is the first one (i.e., 0 ≤ n ≤ 1 ∧ top = 1) then we increment the counter n of detected period, we reset the flag top and update the value of the average duration of detected period tp while we do not update the variable s 2 tp as in order as in order to update the value of the fluctuation of the detected period duration we need that at least two periods have been detected. Finally if the period to be detected is the n-th with n ≥ 2 (i.e., corresponding to guard 2 ≤ n ≤ N ∧ top = 1) we do the same update operations of the previous case but also update s 2 tp . The automata uses variable n to count the number of noisy periods detected along a trajectory, and stops as soon as the N th period is detected (i.e. event bounded measure). The boolean variable top, which is set to true on entering of the high location, allows for detecting the completion of a period (i.e. crossing from mid to low when top is true). Two clock variables, t and t p , maintains respectively the total simulation time as of the beginning of the first detected period (t) and the duration of the last detected period (t p ). Finally variable t p maintain the average duration of all (so far) detected periods while s 2 tp stores the fluctuation (or variability) of duration (i.e. how far the duration of each detected period is distant from its average value computed along a trajectory) of all (so far) detected periods.

Theorem 1 If a trace σ A is noisy periodic w.r.t. amplitude levels L, H ∈ N then it is accepted by automaton A per with parameters L, H, initT ∈ R + and N ∈ N.

Proof See Appendix A.

HASL expressions associated to A per . We define different HASL expressions to be associated to to automaton A per .

-Z 1 ≡ E[last( tp )]: corresponding to the mean value of the period duration for the first N detected periods. -Z 2 ≡ P DF ( tp , s, l, h): corresponding to the PDF of the average period duration over the first N detected periods, where [l, h] represents the considered support of the estimated PDF, and [l, h] is discretized into uniform subintervals of width s -Z 3 ≡ E[last(s 2

tp )]: corresponding to the fluctuation of the period duration.

Expression Z 1 represents the expected value assumed by variable tp , that is, the average duration of the first N periods detected along a trace, at the end of accepted trajectory (i.e., a trajectory that contains N periods). Similarly expression Z 2 evaluates the PDF of the average duration of the first N periods by assuming the interval [l, h] as the support of the PDF and considering that [l, h] is discretised in (hl)/s uniform subintervals of width s. On the other hand Z 3 is concerned with assessing the expected value that variable s 2 tp has at the end of a trace consisting of N noisy periods. By definition (see Table 1) s 2 tp corresponds to the fluctuation of the duration of the detected periods, (i.e., how much the N periods detected along a trace differ from their average duration). Observe that the measured period fluctuation (i.e. Z 3 ) provides us with a useful measure of the irregularity, from the point of view of the period duration, of the observed oscillation.

Measuring the peaks of oscillations

In the previous section we have seen how a characterisation of periodicity for stochastic oscillation can be obtained by considering a given partition, induced by two thresholds L, H, of the domain of the observed population. The drawback of such a characterisation is that, the detected periods depend on the chosen L, H thresholds, and these have to be chosen by the modeller manually, i.e., normally by looking at the shape of a sampled trajectory and then choosing where to "reasonably" set the L and H values before executing the measurements with automaton A per . To improve things here we propose a different approach which is aimed at identifying where the peaks (i.e., the local maxima/minima) of oscillatory traces are located.

Since traces of a DESP consist of discrete increments/decrements of at least one unit, it is up to the observer to establish what should be accounted for as a local maximum (minimum) during such detection process. Intuitively a local max/min of a trace σ A (the projection of σ w.r.t. the observed species A) is a state σ A [i] (i ∈ N) that corresponds to a change of trend in the population of A. This is formally captured by the following definition.

Definition 10 (local maximum/mininimum of a trace) For σ

A the A projection of a trace σ of an n-dimensional DESP D population model, state σ A [i] is a maximum, if σ A [i -1] < σ A [i] > σ A [i + 1], or a minimum, if σ A [i -1] > σ A [i] < σ A [i + 1]. T A m[1]
m [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] m [START_REF] Alur | Model-checking for probabilistic real-time systems[END_REF] m [START_REF] Andrei | Trend-based analysis of a population model of the akap scaffold protein[END_REF] m [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] m [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF] M

[1] M [2] M [3] M [4] M [5] M [6] M [7]
m [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] m [START_REF] Ballarini | Query-based verification of qualitative trends and oscillations in biochemical systems[END_REF] Fig. 5 Example of local maxima/minima of an alternating trajectory σ A .

In the remainder we refer to a that consists of an infinite sequence of local maxima interleaved with an infinite sequence of local minima as an alternating trace (Definition 11).

Definition 11 (alternating trajectory) A trajectory σ of an n-dimensional DESP D population model is said alternating with respect to the i th (1 ≤ i ≤ n) observed species of D, if σ i contains infinitely many local minima (or equivalently local maxima).

For σ A an alternating trace we denote

σ M A = M [1], M [2], . . ., respectively σ m A = m[1]
, m [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF], . . ., the projection of σ A consisting of the local maxima, respectively minima, of σ A . Figure 5 shows the local maxima and minima for an example of alternating trace σ A . In the following we point out two simple properties relating the definition of noisy periodic and alternating trace.

Proposition 1 If σ A is a noisy periodic trace (as of Definition 6) then it is also alternating.

Proposition 1 is trivially true as by definition a noisy period trace visit infinitely often the low and high region of the state space, thus necessarily it contains an infinite sequence local maxima interleaved with local minima.

Proposition 2 If σ A is an alternating trace (as of Definition 11) then it is not necessarily noisy periodic.

Proposition 2 simply points out that, by definition, an alternating trajectory may be diverging (for example if it consist of increasing steps which are always larger than the decreasing ones), in which case clearly it is not noisy periodic.

In the remainder we introduce a HASL based procedure for detecting the local maxima and local minima of alternating traces. However rather considering detection of "simple" local maxima/minima as in Defintion 10, we refer to detection of a generalised notion of local maxima/minima of a trace, that is, maxima and minima which are distanced, at least, by a certain value δ. We formalise this notion in the next definition.

Definition 12 (δ-separated local maxima) Let δ ∈ R + , and σ A the A projection of a trace σ of an n- dimensional DESP D population model. A state σ A [i] is the j-th, j ∈ N >0 , δ-separated local maxi- mum (minimum), denoted M δ [j] (m δ [i]), if it is the largest local maximum (smallest local minimum) whose distance from the preceding δ-separated lo- cal minimum m δ [j -1] (maximum M δ [j -1])
is at least δ and from which smaller (larger) state is reachable which is

> δ distanced from M δ [j] (m δ [i]).
To understand the meaning of δ-separated minima and maxima definition let us consider the example in Figure 6. The nature of the first δ-separated point (either a maximum or a minimum) depends on the initial profile of the trace: if the trace, from its initial state σ A [0], first enters the region σ A [0]+ δ than the first δ-separated point will be a maximum (conversely if it first enters the σ A [0]δ region it will be a minimum). For the trace in Figure 6 the first δ-separated point is the maximum M δ [START_REF]COSMOS home[END_REF] which is preceded by local minima m [START_REF]COSMOS home[END_REF],m [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] and a local maximum M [1] none of which is δdistanced form the initial state σ A [0]. Observe that the detection of M δ [START_REF]COSMOS home[END_REF] as the actual first δ-separated maximum is completed only on entering the M δ [1]δ, which happens on observing the decrement following the local maximum M [START_REF] Alur | Model-checking for probabilistic real-time systems[END_REF]. To better understand that detection of δ-separated minima/maxima may involve a "temporary detection" phase followed by a "detection finalisation" let us continue the unfolding of the trace. Having entered the M δ [START_REF]COSMOS home[END_REF]δ region we encounter the local minimum m [START_REF] Andrei | Trend-based analysis of a population model of the akap scaffold protein[END_REF], which is temporarily detected as the first δ-separated minimum, since in fact it is distanced more than δ from M δ [START_REF]COSMOS home[END_REF]. However m [START_REF] Andrei | Trend-based analysis of a population model of the akap scaffold protein[END_REF] is followed by m [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] a smaller local minimum reached without exiting the m[4] + δ region. Thus m [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] replaces m [START_REF] Andrei | Trend-based analysis of a population model of the akap scaffold protein[END_REF] as temporary first δ-separated minimum. m [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] becomes the actual first δ-separated minimum only on entering the m[5] + δ region (which in this case corresponds with reaching of state M δ [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF]) since no smaller local minima has been detected before entering m[5] + δ.

Definition 13 (δ-separated alternating trace) A trajectory σ of an n-dimensional DESP D population model is said δ-separated alternating with respect to the i th (1 ≤ i ≤ n) observed species of D, if it contains contains infinitely many δ-separated local minima (and equivalently δ-separated local maxima), where δ ∈ R + .

For σ A a δ-separated alternating trace we denote σ

M δ A = M δ [1], M δ [2], . . ., respectively T A m [2] M [1] M [2] m [3] m[2] M [1] m[3] m[4] M [3] M [4] M [5] M [7]
m [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF] m [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] m [START_REF]COSMOS home[END_REF] Fig. 6 Example of δ-separated local maxima (M δ [i]) and minima (m

δ [i]) of an alternating trajectory σ A σ m δ A = m δ [1], m δ [2]
, . . ., the projection of σ A consisting of the δ-separated local maxima, respectively minima, of σ A . Observe that the δ-separated max/min (Figure 6) are a subset of the "simple" max/min (Figure 5). Furthermore the following property holds: Property 1 For δ = 1 the sequence of δ-separated maxima (minima) of an alternating trace σ A coincides with the list of local maxima (minima), that is:

σ M1 A = σ M A and σ m1 A = σ m A .
The detection of the δ-separated local maxima (minima) for a trace σ A can be described in terms of an iterative procedure through which the list of detected max/min are constructed as σ A unfolds. Such a procedure is formally implemented by the LHA A peaks (Figure 7) which we introduce later on. Here, based on the example illustrated in Figure 6, we informally summarise how detection of δ-separated max/min works. The detection requires storing of the most recent (temporary) δ-separated max (min) into a variable named x M (x m ), while once detection of a δ-separated maximum (minimum) is completed the corresponding variable x M (x m ) is copied into a dedicated list, named Lmax, resp. Lmin (see Table 2), which contains the detected points. To understand how detection works let us consider the trace in Figure 6. The first element encountered is the local minimum m [START_REF]COSMOS home[END_REF] which is then stored into x m = m [START_REF]COSMOS home[END_REF]. As the trace further unfolds the subsequent maxima (green points) are ignored as long as their distance from the temporary minimum x m is less than δ, as is the case with M [START_REF]COSMOS home[END_REF]. Similarly any local minimum m[i] (yellow point) that is encountered after that stored in x m is ignored (e.g., m [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF]), unless it is smaller than x m , in which case x m is updated with the newly found smaller minimum. As σ A unfolding proceeds we find the next local max M [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] which is distant more than δ from the temporary minimum x m : this means that x m currently holds an actual δ-distanced minimum hence its value is appended to Lmin and the procedure starts over, in a symmetric fashion, for the detection of the next maximum.

The rational behind the notion of δ-separated max/min is that for locating the actual peaks of a stochastically oscillating trace it is important to be able to distinguish between the minimal peaks corresponding to stochastic noise, the actual peaks of oscillation. With the δ-separated max/min characterisation we provide the modeller with a means to establish an observational perspective: by choosing a specific value for the δ parameters the modeller establishes how big a level of noise he/she wants to ignore when detecting where the oscillation peaks are located.

In the remainder we introduce the LHA A peaks which formally implements the detection of the δseparated peaks of alternating traces.

The automaton A peaks . We introduce an LHA, denoted A peaks (Figure 7), designed for detecting δdistanced local maxima/minima along alternating traces of a given observed species called A. It requires a parameter δ (the chosen noise level) and the partition of the event set E = E +A ∪ E -A ∪ E =A where E +A (respectively E -A , E =A ) is the set of events resulting in an increase (respectively decrease, no effect) of the population of A.

The rationale behind the structure of A peaks is to mimic the cyclic structure of an alternating trace through a loop of four locations, two of which (i.e. Max and Min) are targeted to the detection of local maxima, resp. minima. The simulated trace yields the automaton to loop between Max and Min hence registering the minima/maxima while doing so. The detailed behavior of A peaks is as follows. Processing of a trace starts with a configurable filter of the initial transient (represented as a box in Figure 7) through which a simulated trace is simply let unfolding for a given initT duration.The actual analysis begins in location start from which

start ṫ : 1 noisyDec ṫ : 1 noisyInc ṫ : 1 Min ṫ : 1 ] , ( x > A + ) ^(n M < N 1 ) , { n M + + , L m a x [ x ] + + , S m a x + = x , x : = A } ] , ( x > A ) , { x : = A } E + A , > , ; E =A ,>,; Max ṫ : 1 E +A ,>,{x := A} E =A ,>,; ] , ( x < A ) , { x : = A } E,>,; E,>,;
] ,( x < A

) ,{ x := A } ] ,( x > A + ) ,{ x := A } end ],t = iniT {x := A, n = 1}
A peaks E,>,;

E A ,>,{x := A} E A , > , ; l 0 ṫ : 1 E, >, ;
Skipping the initial transient until iniT time units have passed N n array of frequency of heights of detected maxima (minima)

],(x<A )^(n m = N 1), {Lmin[x]++, Smin+= x, n m ++} ],(x>A+ )^(n M = N 1), {Lmax[x]++, Smax+= x, n M ++} ] ,( x < A ) ^(n m < N 1 ) , { n m + + , L m i n [x ] + + , S m i n + = x , x := A }
Table 2 The data variables of automaton A peaks of Figure 7 for locating the peaks of a noisy oscillatory traces we move to either Max or Min depending whether we initially observe an increase (i.e. x < Aδ) or a decrease (i.e. x > A + δ) of the population of the observed species A beyond the chosen level of noise δ. Once within the Max→noisyDec →Min→ noisyInc loop the detection of local maxima and minima begins. Location Max (Min) is entered from noisyInc (noisyDec) each time a sufficiently large (w.r.t. δ) increment (decrement) of A is observed. On entering Max (Min), we are sure that the current value of A has moved up (down) of at least δ from the last value stored in x while in Min (Max), hence that value (x) is an actual local minimum (maximum) thus we add it up to Smin (Smax), then we increment the frequency counter corresponding to the level of the detected mini-mum Lmin[x] (maximum Lmax[x]) 10 before storing the new value of A in x and finally increase n M (n m ) the counter of detected maxima (minima). Once in Max (Min) we stay there as long as we observe the occurrence of reactions which do not decrease (increase) the value of A, hence either reactions of E +A (E -A ) or of E =A . While in Max (Min) if we observe a reaction of E +A (E -A ), then we store the new increased (decreased) value of A in x, which becomes the new potential next local maximum (minimum). On the other hand while in Max (Min) if a "decreasing" ("increasing") reaction E -A (E +A ) occurs we move to noisyDec (noisyInc) from which we can either move back to Max (Min), if we observe a new increase (decrease) that makes the population of A overpass x (x overpass A), or eventually entering Min (Max) as soon as the observed decrease (increase) goes beyond the chosen δ (see above). For the automaton A peaks depicted in Figure 7, the analysis of the simulated trace ends, by entering the end location either from noisyDec or noisyInc, as soon as N maxima (or minima, depending on whether the first observed peak was a maximum or a minimum) have been detected. Notice that A peaks can straightforwardly be adapted to different ending conditions. The data variables of A peaks are summarised in Table 2.

HASL expressions associated with A peaks . We define different HASL expressions to be associated with automaton A peaks .

-Z max ≡ E[last(Smax)/n M ]: corresponding to the expected value of the average height of the maximal peaks for the first N detected maxima. 10 with a slight abuse of notation we refer to Lmin[] and Lmax[] as arrays whereas in reality within COS-MOS/HASL they correspond to a set of variables Lmin i , Lmax j , each of which is associated to a given level of the observed population, thus Lmin 1 counts the frequency of observed minimum at value 1, Lmin 2 the observed minima at value 2 and so on. The number of required Lmin i , Lmax j variables, which is potentially infinite, can be actually bounded without loss of precision to a sufficiently large value Lminm (resp. Lmaxm) which must be established manually beforehand, for example by observing few previously generated traces.

-Z min ≡ E[last(Smin)/n m ]: same as Z max but for minima. -Z P DF max ≡ E(last(Lmax)/n M ): enabling to compute the PDF of the height (along a path) of the maximal peaks -Z P DF min ≡ E(last(Lmin)/n m ): enabling to compute the PDF of the height (along a path) of the maximal peaks Expression Z max (Z min ) represents the average value of the detected δ-separated maxima (minima). This is obtained by considering the sum of all detected δ-separated local maxima (minima), which is stored in Smax (Smin) and dividing it by the number of detected maxima n m (n m ). Expression Z P DF max (Z P DF min ) allows to estimate the PDF of the height of the detected δ-separated local maxima (minima). This is achieved by dividing the frequency counters of each detected maximal (minimal) peak's height, whose values are stored in array Lmax (Lmin), by n M (n m ), the number of detected maxima (minima).

Theorem 2 If σ A is a δ-separated (with δ ∈ R + ) alternating trace then it is accepted by automaton A peaks with parameters δ, initT ∈ R + and N ∈ N.

Proof See Appendix A.

Case study

To demonstrate the above described procedure we consider a popular example of oscillator, the socalled circadian clock. Circadian clocks are biological mechanisms responsible for keeping track of daily cycles of light and darkness. Here we focus on a model of the biochemical network ( [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF]) which is believed to be at the basis of the control of circadian clocks. The network (Figure 8 by gene D R . M A and M R are then translated into proteins A, respectively R .

R 1 :A + D A γ A -→ D A R 9 : M A β A -→ M A + A R 2 :D A θa -→ A + D A R 10 : M R β R -→ M R + R R 3 :A + D R γ R -→ D R R 11 : A + R γ C -→ C R 4 :D R θ R -→ D R + A R 12 : C δ A -→ R R 5 :D A α A -→ M A + D A R 13 : A δ A -→ ∅ R 6 :D A α A -→ M A + D A R 14 : R δ R -→ ∅ R 7 :D R α R -→ M R + D R R 15 : M A δ M A -→ ∅ R 8 :D R α R -→ M R + D R R 16 : M R δ M R -→ ∅ (2)
Protein A acts as an activator for both genes by attaching to promoter region of D A and D R (i.e. when A is attached to a gene the mRNA transcription increases). Species D A and D R represent the state of gene D A , respectively D B , when an activator molecule (A) is attached to their promoter. Note that gene D R acts as a repressor of D A since when A bounds to its promoter D R sequesters the activator A and, as a result, the transcription of D A slows down. The repressing role of D R is further due to the fact that the expressed protein R inactivates the activator A by binding to it and forming the complex C. Finally the model in Figure 8 accounts for degradation of all species: thus the mR-NAs M A and M R , as well as the expressed proteins A and B degrades with given rates (see Table 3). Notice that protein A degrades also when attached to R (i.e. when in complex C), and, as a consequence, C turns into R at a rate equivalent to the degradation rate of A. The model of Figure 8 corresponds to the system of chemical equations (2), whose (continuous) kinetic rates (taken from [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF]) are given in Table 3.

Stochastic model Equations ( 2) can give rise to either a system of ODEs or to a stochastic process.

Here we focus on the discrete-stochastic semantics: Figure 9 shows the GSPN encoding of equations (2) developed with COSMOS. The configuration of the GSPN (i.e. the stochastic process) requires setting the initial population and the rates of each transition (i.e. reaction). For the initial population, following [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF], we observe that the model comprises one gene D A and one D R , which can either be in free-state (no activator A is attached to the promoter) or in activator-bound state, i.e. D A , respectively D R . As a consequence the population of species D A and D R is bounded by the following invariant constraints: D A + D A = 1 and D R + D R = 1 (in fact places DA, DA and DR, DR of net in Figure 9 are the only places covered by P-invarriants). The remaining species are initially supposed to be "empty", hence they are initialised to 0. Concerning the transition rates, for simplicity we assume a unitary volume of the system under consideration, hence all continuous rates in Table 3 can be used straightforwardly as rates of the corresponding discrete-stochastic reactions. In this case we assume all reactions following a negative exponential law.

The oscillatory dynamics of the GSPN model of Figure 9 can be observed by plotting of a simulated trajectory (Figure 10). Observe that the frequency of oscillations varies considerably with the degradation rate of the repressor (R) protein: a faster degradation of R (right), intuitively, results in a higher frequency of oscillations. In the remainder we formally assess the oscillatory characteristics (i.e. the period and the peaks of oscillations) of the circadian clock model by application of the previously described approach, i.e. by analysing the stochastic process deriving from synchronisation of the circadian clock GSPN model with the A period and A peaks automata.

α A 50 h -1 α R 0.01 h -1 δ A 1 h -1 δ R 0.2 h -1 α A 500 h -1 α R 50 h -1 γ A = γ R 1 mol -1 h -1 γ C 2 mol -1 h -1 β A 50 h -1 β R 5 h -1 θ A 50 h -1 θ R 100 h -1 δ M A 10 h -1 δ M R 0.5 h -1
Table 3 reactions' rates for the circadian oscillator

R 15 M A M A R 6 ↵ A R 9 A D A R 1 A R 2 ✓ A R 4 ✓ R R 3 R D R R 8 ↵ R M R R 16 M R R 11 C A R R 1 R C R 5 ↵ 0 A D 0 A D 0 R R 7 ↵ 0 R R 12 A R 14 R 1
Fig. 9 GSPN encoding of the system (2) of chemical equations corresponding to the circadian-clock. Measuring the period of the circadian clock. We performed a number of experiments aimed at assessing the effect that the degradation rate of the repressor protein (δ R ) has on the period of the circadian oscillator. Figure 11 (right) shows three plots representing the PDF of the period (obtained through the HASL formula (A period , P DF (Last(t)/N )) for three values of δ R . With δ R = 0.2 (i.e. the original value as given in [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF]) the PDF is centred at t = 24.9, i.e. slightly more of the standard 24 hours period expected for a circadian clock. On the other hand speeding up the repressor degradation of 10 times (i.e. δ R = 2) yields a slightly more than halved oscillation period (i.e. PDF cen-tred at T = 10.8). Finally slowing down the degradation rate of a half (i.e. δ R = 0.1) yields a less than doubled oscillation period (i.e. PDF centred at T = 40.7). where periods have been detected using L = 1 and H = 1000 as partition thresholds, and target estimates have been computed with confidence level 99 and confidence-interval width of 0.01. Furthermore the PDF plots in Figure 11 (right) have been computed using a discretisation of the period support interval [0, 50] into subintervals of width 0.1.

Measuring the peaks of oscillations of the circadian clock. We performed a number of experiments aimed at assessing the effect that the degradation rate of the repressor protein (δ R ) has on the peaks of oscillation for both protein A and R. Figure 12 shows plots for the mean value of the minimal and maximal peaks of oscillations for both A and R. Results indicate that while the degradation rate δ R has no effect on the oscillation peaks of A (both maximal and minimal peaks of A are constant independently of δ R ), it affects the maximal peaks (only) of R. Specifically the mean value of R's maximal peaks decreases with the increasing of δ R (while the minimal peaks of R are constantly at 0), notice that this is in agreement with what indicated by the single trajectories depicted in Figure 10. Notice that Figure 12 contains also plot for the absolute maximum of population of A and R measured along the sampled trajectories through trivial HASL expressions Z ≡ AV G(max(x)) (where x is a variable used to record the population of the observed species along a synchronising path). All plots in Figure 12 result from sampling of finite trajectories containing of N = 100 maximal peaks and using a noise parameter δ = 10%AV G(max(x)), meaning that for evaluating the mean value of maximal peaks we discarded all critical points distanced one another less than 10% of the absolute maximum of the observed species. Finally, again points of every plot in Figure 12 have been computed with confidence level 99 and confidence-interval width of 0.01.

On the initial transient. To assess the effect that the initial transient of the circadian clock model have on the period and peaks estimates we repeated all of the above discussed experiments with different values of the initT parameter (e.g. initT ∈ {10, 50, 100, 500, 1000, . . .}) which determines the starting measuring point for A period and of A peaks . The outcomes of repeated experiments turned out to be independent of the chosen initT value, indicating that circadian clock reaches its steady state very quickly.

Related work and discussion

The HASL based methodology presented in this paper is by no means the only approach aimed at the analysis of discrete-state stochastic oscillators. In the following we provide a brief (non exhaustive) overview of similar approaches.

Mathematical approaches

The analysis of periodic signals can be achieved through well established signal processing techniques such as, for example, Fast Fourier Transform (FFT) and autocorrelation. Both methods estimate the dominant frequency of a periodic signal given in terms of a sequence of (real-valued) points. In the context of stochastic modelling both FFT and autocorrelation analysis is performed over trajectories generated by a stochastic simulator. In order to increase the accuracy of the estimates usually frequency estimation is then replicated over N trajectories, the final result being given as the average of the frequency estimate of each trajectory (see e.g. [START_REF] David | Runtime verification of biological systems[END_REF][START_REF] Ihekwaba | Communicating oscillatory networks: frequency domain analysis[END_REF]). The main appeal of signal processing techniques is due to their simplicity. However, in the context of statistical model checking, adding an (automatic) control on the accuracy of the resulting estimate would require their integration within a confidence interval estimation procedure, something which at best of our knowledge has not yet been done. From an expressiveness point of view it is worth remarking that FFT and autocorrelation are limited to estimating the (mean value) of the frequency of an oscillator but provides no support for assessing other aspects of oscillator such as the location of the oscillation peaks and the regularity (i.e. the fluctuation) of the period. Finally another interesting contribution belonging to the field of mathematical approaches is presented in [START_REF] Júlvez | Evaluation of sustained stochastic oscillations by means of a system of differential equations[END_REF], where the relationship between stochastic oscillators and their continuous-deterministic counterpart is analysed.

Model checking based approaches. Analysis of oscillators through stochastic model checking techniques has been considered in several works. Application of CSL [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] to the characterisation CTMC biochemical oscillators has been considered, with limited success in [START_REF] Ballarini | Query-based verification of qualitative trends and oscillations in biochemical systems[END_REF], and more comprehensively in [START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in markovian population models[END_REF][START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in markov population models[END_REF]. In [START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in markovian population models[END_REF] Spieler demonstrated that deciding whether a given CTMC model oscillates sustainably boils down to a steady-state analysis problem where the allegedly oscillating CTMC is coupled with a period detector automata (through manual hard-wiring). In this case the probability that the period of oscillation has a certain value is computed through dedicated CSL steady-state formulae and has been demonstrated through examples on the PRISM model-checker.

In a recent work measuring of oscillations has been considered with other statistical model checking tools (UPPAAL-SMC and PLASMA) by application of the MITL logic [START_REF] David | Runtime verification of biological systems[END_REF]. In this case the analysis of period duration is achieved by detection of a single period of oscillation through nested time-bounded Until formulae.

Conclusion

We have presented a methodology for the formal analysis of stochastic models exhibiting an oscillatory behaviour. Such a methodology relies on the application of the HASL formalism, a statistical model checking framework suitable for expressing sophisticated performance measures. We have shown how by means of HASL one can define specific LHA automata targeted to the analysis of particular aspects of the dynamics of oscillatory trajectories, such as the detection of the period and of the peaks of a stochastic oscillator. For the period we have introduced a template LHA, denoted A per , which allows us to estimate the PDF, the mean value as well as the the fluctuation of the period duration. Concerning the peaks we have introduced a template LHA, A peaks , which allows us to measure the mean value of the maximal/minimal peaks of oscillation. We have demonstrated the effectiveness of the methodology by studying a well established model of biological oscillator, namely the circadian A Appendix Theorem 1. If a trace σ A is noisy periodic w.r.t. amplitude levels L, H ∈ N then it is accepted by automaton Aper with parameters L, H, initT ∈ R + and N ∈ N Proof By hypothesis σ A (the projection, w.r.t an observed species A, of a trace σ an n-dimensional DESP D) is noisy periodic w.r.t. the partition of species A's domain into regions low = (-∞, L), mid = [L, H) and high = [H, ∞).

The initial state of the synchronised process D × Aper will be (σ A [0], l 0 , V al 0 ) with V al 0 being the initial valuation with V al 0 (x) = 0 for all variables ∀x ∈ X of Aper.

To demonstrate that σ A is accepted by Aper we need to show that starting from the the initials state (σ A [0], l 0 , V al 0 ) a final state of D × Aper, i.e., a state such that the current location is the accepting location end of Aper, is reached. For this we proceed by induction w.r.t. the number of subsequent sojourns in the low and high regions. In the remainder we use the notation (s, l, V al) * -→ (s , l , V al ) to indicate that state (s , l , V al ) of process D × Aper is reachable from (s, l, V al). We split the demonstration in parts corresponding to the traversal of Aper locations resulting from synchronisation with trace σ A :

[init ] Let l 0 ∈ N be the index of the first state of σ A that belongs to low and that follows σ A @initT (where initT is the parameter of Aper), that is:

l 0 = min{i ∈ N | i > i initT ∧ σ[i] ∈ low}, with i initT being the index of the state σ A is in at time initT (observe that since σ A is assumed noisy periodic then min{i ∈ N | i > i initT ∧ σ[i] ∈ low} is guaranteed to exist). Thus because of the structure of Aper, (σ A [0], l 0 , V al 0 ) * -→ (σ[l 0 ], low, V al l0 ), with V al l0 (x) = 0, for x = n and V al l0 (n) = -1.
[low→mid→high ] Since σ A is noisy periodic then ∃m 1 , h 1 ∈ N : [high→mid→low ] similarly since σ A is noisy periodic then ∃l 1 , m 1b ∈ N :

h 1 > m 1 > l 0 such that σ A [m 1 ] ∈ mid, σ A [h 1 ] ∈ high hence,
l 1 > m 1b > h 1 such that σ A [m 1b ] ∈ mid, σ A [l 1 ] ∈ low. hence (σ[h 1 ], high, V al h1 ) * -→ (σ[m 1b ], mid, V alm 1b ) * -→ (σ[l 1 ]
, low, V al l1 ) with V alm 1b (top) = 1, V alm 1b (n) = -1 hence V al l1 (top) = 0, V al l1 (n) = 0, V al l1 (t) = 0, since, because of V alm 1b , location low is entered through edge mid E,(n=-1∧top=1),{n++,top:=0,t:=0} --------------------------→low induction since σ A is noisy periodic then the high, low regions are entered infinitely often, hence the [low→mid→high] and [high→mid→low] steps of the proof hold for each successive iteration. This means that if l i is the index corresponding to the ith time that σ A enters the low region after having sojourned in the high region then because of the periodicity of σ A ∃m i+1 , h i+1 , m(i + 1)b ∈ N : Theorem 2. If σ A is a δ-separated (with δ ∈ R + ) alternating trace then it is accepted by automaton A peaks with parameters δ, initT ∈ R + and N ∈ N.

l i+1 > m (i+1)b > h 1 > m 1 > l i such that σ A [m ib ], σ A [m (i+1)b ] ∈ mid, σ A [l i ], σ A [l i+1 ] ∈ low, σ A [h (i+1) ] ∈ high, hence (
Proof By hypothesis σ A is δ-separated alternating, which means it contains an infinity of δ-separated local maxima interleaved with δ-separated minima. We recall the basic notation: if σ A is alternating then Let (σ A [0], l 0 , V al 0 ) be the initial state of the process corresponding to the synchronization of σ A with A peaks , with V al 0 being the initial valuation with V al 0 (x) = 0 for all variables ∀x ∈ X of A peaks . We need to show that a final state (σ A [n], end, V aln) is reached when σ A is synchronised with A peaks . For this we proceed by induction w.r.t. the number of subsequent δ-separated max and min encountered during the synchronization of δ A . We split the demonstration in parts corresponding to the traversal of A peaks locations resulting from synchronisation with trace σ A .

[l 0 → start] Because of the structure of A peaks , we have that (σ A [0], l 0 , V al 0 ) * -→ (σ A @initT, start, V alstart), with V alstart(n) = -1, V alstart(x) = A and V alstart(t) = initT .

[start→Max or start→Min] Let i 0 ∈ N be the index of the first state of σ A that follows σ A @initT and that is distanced more than δ from σ A @initT , i.e., i 0 = min{i ∈ N | σ A @initT -σ A [i 0 ] > δ}, (observe that since σ A is assumed to be δ-separated alternating then min{i ∈ N | σ A @initT -σ A [i 0 ] > δ} is guaranteed to exist). Furthermore let i initT be the index of state σ A @initT , i.e. σ A [i initT ] = σ A @initT , and ∆T = n|i initT +n<i0 i=0 time(σ A , i initT + i), be the time distance from initT and the instant when state σ A [i 0 ] is entered. Thus because of the structure of A peaks , either (σ A @initT, start, V alstart) * -→ (σ[i 0 ], Max, V al i0 ), if (σ A @initT-σ A [i 0 ]) < 0, or (σ A @initT, start, V alstart) * -→ (σ[i 0 ], Min, V al i0 ), if (σ A @initT -σ A [i 0 ]) > 0, with V al i0 (x) = A, V al i0 (t) = V alstart(t) + ∆T .

[Max→NoisyDec→Min ] Let us assume that from state (σ A @initT, start, V alstart) we reached (σ[i 0 ], Max, V al i0 ), which means that from σ A @initT the trace σ A reaches a state σ A [i 0 ] > σ A @initT + δ without ever descending below σ@initT -δ, i.e. ∀j < i 0 , σ A [j] ≥ σ A @initTδ. [Min→NoisyInc→Max ] On the contrary if we assume that from state (σ A @initT, start, V alstart) we reached (σ[i 0 ], Min, V al i0 ), which means that from σ A @initT the trace σ A reaches a state σ A [i 0 ] < σ A @initT -δ without ever entering the σ@initT +δ region, i.e. ∀j < i 0 , σ A [j] ≤ σ A @initT +δ. Since σ A is δ-separated alternating and since we assumed to have reached configuration (σ[i 0 ], Min, V al i0 ) then the first δ-separated point is a minimum and in particular is a state σ A [h 0 ] = m δ [1], i 0 < h 0 , such that σ A [h 0 ] is the smallest state reachable from σ A [i 0 ] (i.e., σ A [j ] ≥ σ A [h 0 ], ∀j , i 0 < j < h 0 ) and from which a δ-distanced larger state

[k 0 ], i 0 < H 0 < k 0 is reachable with k 0 = inf {k | σ A [k] < σ A [H 0 ] - δ ∧ σ A [k ] ≤ σ A [H 0 ], ∀k , H 0 < k < k}.
σ A [k 0 ], i 0 < h 0 < k 0 is reachable with k 0 = inf {k | σ A [k] > σ A [h 0 ] -δ ∧ σ A [k ] ≥ σ A [h 0 ], ∀k , h 0 < k < k}.
Thus, because of the structure of A peaks , it follows that (σ[i 0 ], Min, V al i0 induction since σ A is δ-separated alternating then it contains and infinity of δ-separated maxima and minima. hence the [Max→NoisyDec→Min] and [Min→NoisyInc→Max] steps of the proof hold for each successive iteration. This means that if k i , K i ∈ N are the index corresponding to the state of σ A determining the detection of the i-th δ-separated maximum, resp. minimum, then ------------------------

Definition 3 (

 3 DESP Population Model) A DESP model for n ∈ N population types is a DESP D = S, π 0 , E, Ind, enabled, delay, choice, target with S ⊆ N n . Definition 4 (DESP Observed Species) For D a DESP population model with n species we define D i the observed i th process, with 1 ≤ i ≤ n, as the process resulting from D by observing only the i th component of each state of D. Thus each s

EFig. 1

 1 Fig. 1 Synchronisation between a DESP toy model (left) representing basic steps of gene expression and an LHA (right) which selects paths containing N occurrences of the transcription event

-φ 1 ≡

 1 (A, E[last(t)]): represents the average time for completing N transcriptions φ 2 ≡ (A, E[max(a)]): represents the maximum population reached by protein A within the first N transcriptions φ 3 ≡ (A, P DF (last(t), 0.1, 0, 10)): represents the PDF of the delay for completing N transcriptions computed over the interval [0, 10] with a discretisation step of 0.1

Fig. 2

 2 Fig. 2 Deterministic versus stochastic (noisy) oscillations

EFig. 4

 4 Fig.4Aper: an LHA for selecting noisy periodic traces (with respect to an observed species A) related to partition low = (-∞, L], mid = (L, H) and high = +∞).

Fig. 7 A

 7 Fig.7A peaks : an LHA for detecting local maxima/minima (for observed species A) of noisy periodic traces where local maxima/minima are detected with respect to a chosen level of noise δ.Data variables name domain update definition descriptiont

Fig. 8

 8 Fig. 8 Circadian Clock oscillator network: gene D A expresses activator protein A through transcription of mRNA M A , gene D R expresses the repressor protein R through transcription of mRNA M R .

Fig. 10

 10 Fig. 10 Single trajectory showing the oscillatory character of activator A and repressor B dynamics with normal repressor's degradation rate δ R = 0.2 (left) and with 10× speed-up, i.e. δ R = 2 (right).

Figure 11 (Fig. 11 Fig. 12

 111112 Figure 11 (left) shows plots for the period mean value (red plot) and the period fluctuation (blue plot, as described in Section 3.1) in function of the degradation rate δ R . They indicate that slowing down the degradation of the repressor yields, on one hand, to a lower the frequency of oscillations, and on the other, augmenting the irregularity of the periods (i.e. augmenting the period's fluctuations). All plots in Figure11result from sampling of finite trajectories consisting of N = 100 periods,

  because of the structure of Aper, it follows that (σ[l 0 ], low, V al l0 ) * -→ (σ[m 1 ], mid, V alm 1 ) * -→ (σ[h 1 ], high, V al h1 ) with V al h1 (top) = 1,because of the update {top := 1} of the arc leading to location high.

  σ[l i ], low, V al l i ) * -→ (σ[l i+1 ], low, V al l i+1 ), with V al l i+1 (n) = V al l i (n) + 1. termination, [low→end ] By induction we have seen that ∀i ∈ N, (σ[l i ], low, V al l i ) * -→ (σ[l i+1 ], low, V al l i+1 ). Thus on the (N -1)-th iteration (σ[l N -1 ], low, V al l N -1 ) * -→ (σ[l N ], low, V al l N ) with V al l N (n) = N which enables low ,(n=N ),∅ -------→end, hence (σ[l N ], low, V al l N ) * -→ (σ[l N ], end, V al l N ) and σ A is accepted.

σ

  A = * , m[1], * , M [1], * , m[2], * , . . . where m[i] (res[. M [i]) are the local minima (resp. maxima) and * denotes an arbitrary long sequence of states.

A

  = M[START_REF]COSMOS home[END_REF], M[START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF], . . ., respectively σ m A = m[START_REF]COSMOS home[END_REF], m[START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF], . . ., is the projection of σ A consisting of the local maxima, respectively minima, of σ A . If σ A is δ-separated alternating then:σ A = σ 1 A , m δ [1], σ 2 A , M δ [1], . . . where each σ i A (i ∈ N) is a finite alternating sequence of non-δseparated local minima/maxima and m δ [i] (resp. M δ [i]) are the δ-separated local minima (resp. maxima) Forδ ∈ R + we denote σ M δ A = M δ [1], M δ [2], . . ., respectively σ m δ A = m δ [1], m δ [2], . . ., the projection of σ A consisting of the δ-separated local maxima, respectively minima, of σ A .

  Thus, because of the structure of A peaks , it follows that (σ[i 0 ], Max, V al i0 ) * -→ (σ[k 0 ], Min, V al k0 ), where V al k0 (Lmax[σ[H 0 ]]) = V al i0 (Lmax[σ[H 0 ]]) + 1 V al k0 (n M ) = V al i0 (n M ) + 1, V al k0 (x) = σ A [k 0 ], V al k0 (Smax) = V al i0 (Smax) + 1, V al k0 (t) = V al i0 (t) + ∆T , with ∆T = k0j=i0 time(σ A , j).

  ) * -→ (σ[k 0 ], Max, V al k0 ), where V al k0 (Lmin[σ[h 0 ]]) = V al i0 (Lmin[σ[h 0 ]]) + 1, V al k0 (nm) = V al i0 (nm) + 1, V al k0 (x) = σ A [k 0 ], V al k0 (Smin) = V al i0 (Smin)+1, V al k0 (t) = V al i0 (t)+∆T , with ∆T = k0 j=i0 time(σ A , j).

  (n M ) = V al k i (n M + 1). Thus on he (N M -1)-th iteration we have V al k n M -1 = N M -1, hence on detection of the N M -th δ-separated maximum the LHA edge NoisyDec ,(x>A+δ)∧(n M =N M -1)

  → {Lmax[x]++,Smax+=x,n M ++} end is enabled hence (σ[k n M -1 ], Max, V al k n M -1 ) * -→ (σ[Kn M ], end, V al Kn M ), and σ A is accepted. termination, [Min→NoisyInc→end ] specular to [termination, [Max→NoisyDec→end] ] .

  Algorithm 1 Evolution of a DESP initial configuration: (s , τ , sched) while Enabled(s) = ∅ do E

min = min(∪ e∈enabled(s) sched(e)) enext = choice(s, E min , τ ) δm = sched(enext) s = target(s, enext, δm) τ = τ + δm sched(e) = +∞ (∀e ∈ enabled(s )) sched(e) = sample(delay(s , e)) (∀e ∈ enabled(s ))

  Since σ A is δ-separated alternating and since we assumed to have reached configuration (σ[i 0 ], Max, V al i0 ) then the first δ-separated point is a maximum and in particular is a stateσ A [H 0 ] = M δ [1], i 0 < H 0 , such that σ A [H 0 ] is the largest state reachable from σ A [i 0 ] (i.e., σ A [j ] ≤ σ A [H 0 ], ∀j , i 0 < j < H 0 )and from which a δ-distanced smaller state σ A

  ∀i ≥ 0 either (σ[K i ], Min, V al K i ) * -→ (σ[k i ], Max, V al k i ) * -→ (σ[K i+1 ], Min, V al K i+1 ), if (σ A @initT, start, V alstart) * -→ (σ[i 0 ], Max, V al i0 ), or (σ[k i ], Max, V al k i ) * -→ (σ[K i ], Min, V al K i ) * -→ (σ[k i+1 ], Max, V al k i+1 ), if (σ A @initT, start, V alstart) * -→ (σ[i 0 ], Min, V al i0 ), given that V al K i (nm) < Nm -1 and V al k i (n M ) < N M -1.termination, [Max→NoisyDec→end ] By induction we have seen that ∀i ∈ N, given thatV al k i (n M ) < N M -1, we have (σ[k i ], Max, V al k i ) * -→ (σ[K i ], Min, V al K i ) where V al K i

hence, in essence, a DESP corresponds to a generalised semi-Markov processes[14,

[START_REF] Alur | Model-checking for probabilistic real-time systems[END_REF]

GSPN with timed transitions associated to generic probability distributions, that is, not necessarily Negative Exponential as with the standard GSPN definition[START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF].

Notice that because of the "initial-nondeterminism" of LHA there can be at most one initial state for the product process.

notice that because of the determinism constraints of the LHA edges (conditions c2 and c3) at most only one autonomous or synchronised edge can ever be enabled in any location of the LHA.

Otherwise PDF/CDF measures can be encoded explicitly in an LHA but such encoding would usually result in a rather complex LHA.

COSMOS is an acronym of the french sentence "Concept et Outils Statistiques pour le MOdèles Stochastiques" whose english translation would sound like: "Tools and Concepts for Statistical analysis of stochastic models".

Note that the duration of the k-th period of a trace is, in turn, dependent on the the random variables x j↓ , x j↑ corresponding to the entering of the low, high regions.

Although the LHA in Figure4is designed so that periods detection starts from low it can be easily adapted so that the identification starts from any location.