Paolo Ballarini

Marie Duflot

Applications of an expressive statistical model checking approach to the analysis of genetic circuits

Keywords: model checking, stochastic methods, qualitative/quantitative analysis of biochemical systems, delayed stochastic dynamics

Stochastic temporal logics have demonstrated their efficiency in the analysis of discrete-state stochastic models. In this paper we consider the application of a recently introduced formalism, namely the Hybrid Automata Stochastic Language (HASL), to the analysis of biological models of genetic circuits. In particular we demonstrate the potential of HASL by focusing on two aspects: first the analysis of a genetic oscillator and then the analysis of gene expression. With respect to oscillations, we formalize a number of HASL based measures which we apply on a realistic model of a three-gene repressilator. With respect to gene expression, we consider a model with delayed stochastic dynamics, a class of systems whose dynamics includes both Markovian and non-Markovian events, and we identify a number of relevant and sophisticated measures. To assess the HASL defined measures we employ the COSMOS tool, a statistical model checker designed for HASL model checking.

Introduction

Biological systems are regulated by complex information processing mechanisms which are at the basis of their survival and adaptation to environmental changes. Despite the continuous advancements in experimental methods, many of those mechanisms are still not well understood. The end goal of computational systems biology [START_REF] Kitano | Foundations of Systems Biology[END_REF] is thus to develop formal methods for rigorously representing and effectively analyzing biological systems. Understanding what cells actually compute, how they perform computations, and eventually how such computations can be modified/engineered, are essential tasks which computational modeling aims at. In this context, the ability to "interrogate" a model by asking relevant "questions", referred to as model checking, is critical. Model checking approaches have proved to be an effective means for analyzing biological systems, both in the framework of non-probabilistic models [START_REF] Fages | Modelling and querying interaction networks in the biochemical abstract machine BIOCHAM[END_REF][START_REF] Bernot | Application of formal methods to biological regulatory networks: extending Thomas' asynchronous logical approach with temporal logic[END_REF] and in that of stochastic models [START_REF] Lakin | Design and analysis of DNA strand displacement devices using probabilistic model checking[END_REF][START_REF] Heath | Probabilistic model checking of complex biological pathways[END_REF][START_REF] Barbuti | Probabilistic model checking of biological systems with uncertain kinetic rates[END_REF]. In the realm of stochastic modeling, such questions have a quantitative nature as they most often take into account timing as well as likelihood of events. In this paper we consider the application of a recently introduced stochastic formalism, namely HASL, to the verification of stochastic models of biological mechanisms. HASL uses Linear Hybrid Automata (LHA) as a means for characterizing relevant dynamics of a model, and, as such, it allows one to conceive and evaluate sophisticated measures.

Our contribution. We consider the application of HASL to the verification of biological systems. Our contribution is twofold. In the first part we illustrate the application of HASL to a relevant aspect of many biological mechanisms, namely the analysis of oscillatory trends in stochastic models of biological systems. In so doing we extend the very preliminary results presented in [START_REF] Ballarini | Expressive statistical model checking of genetic networks with delayed stochastic dynamics[END_REF] (where we showed how HASL can be used to count the average number of oscillation periods on a simple example of a 3-species oscillator) by introducing an elaborate approach based on two different oscillation related behaviors (i.e. noisy periodicity and noisy alternance). We then illustrate the proposed approach on a stochastic oscillator case-study: the so called repressilator.

In the second part, we demonstrate the effectiveness of HASL verification by developing a comprehensive case study of gene expression, a relevant biological mechanism represented by means of non-Markovian models, that cannot be analyzed by means of classical (numerical) stochastic model checking.

Paper organisation. In Section 2 we present the background material the paper is based upon, i.e. how systems can be modeled and analyzed using HASL formalism. Section 3 is devoted to the application of HASL to the analysis of oscillations while Section 4 illustrate the HASL-based analysis of a model of gene-expression. Finally a conclusion is given in Section 5.

Statistical model checking with HASL

Model checking [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF] is a well established methodology introduced in the early 80s for the automatic verification of discrete-state models against temporal logic properties. The original approach was then extended to probabilistic models, and in particular to continuous-time Markov Chains (CTMCs) [START_REF] Stewart | Introduction to the numerical solution of Markov chains[END_REF], resulting in the introduction of the Continuous Stochastic Logic (CSL) [START_REF] Aziz | Model-checking CTMCs[END_REF][START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] and subsequent extensions. The basic idea of stochastic model checking is that, given a stochastic model and a property (expressed in terms of a stochastic logic formula), algorithms (either based on numerical methods or statistical ones) are applied to approximate the probability with which the property is fulfilled by the model.

The Hybrid Automata Stochastic Language (HASL) [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] is a recently introduced formalism widening the family of model checking approaches for stochastic models. Its main characteristics are as follows: first it addresses a broad class of stochastic processes which includes but, unlike most stochastic logics, is not limited to CTMCs. Second the HASL formalism turns out to be a powerful specification language through which temporal reasoning is naturally blended with elaborate reward-based analysis. In that respect HASL unifies the expressiveness of CSL [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF] and its action-based [START_REF] Baier | Model checking action-and state-labelled Markov chains[END_REF], timed-automata [START_REF] Donatelli | Model checking timed and stochastic properties with CSL T A[END_REF][START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF] and reward-based [START_REF] Haverkort | Model checking performability properties[END_REF] extensions, in a single powerful formalism. Third HASL model checking belongs to the family of statistical model checking approaches (i.e. those that employ stochastic simulation as a means to obtain estimates of relevant properties of the considered model). More specifically HASL statistical model checking employs confidence-interval methods to estimate the expected value of random variables which may represent either a measure of probability or a generic real-valued measure.

In the following we recall the basics of the HASL formalism. First we present Discrete Event Stochastic Process (DESP), the expressive class of stochastic models to which HASL formalism can be applied. Then we describe an extended version of generalized stochastic Petri nets (GSPN), called non Markovian stochastic Petri nets (NMSPN), which is the modelling formalism used within HASL (note that the semantics of an NMSPN is a DESP). We then introduce the notion of HASL specifications, which consist of an automaton for selecting relevant executions of the DESP and obtaining information about them, and a target measure to be evaluated. We finally briefly describe the statistical model checking scheme, upon which the COSMOS tool [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] (i.e. the model checker supporting the HASL approach) is based. For a comprehensive and more formal treatment of HASL we refer the reader to [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF].

Discrete Event Stochastic Process

A DESP is a stochastic process consisting of a (possibly infinite) set S of states and whose dynamic is triggered by a (finite) set E of discrete events. As mentioned beforehand, no restrictions are considered on the nature of the delay distribution associated with events, thus any distribution with non-negative support may be considered. This model is very similar to Generalized Semi Markov Processes [START_REF] Glynn | A gsmp formalism for discrete event systems[END_REF]. The only difference, set aside the way to present the functions, being the introduction of "indicators" needed when synchronizing the DESP with an LHA. The automaton will indeed not have access to the actual state of the DESP but will get partial information through the values of chosen indicators. Definition 2.1. A DESP is a tuple D = S, π 0 , E, Ind, enabled, target, delay, choice where:

• S is a discrete set of states,

• π 0 ∈ dist(S) is the initial distribution on states,

• E is a set of events,

• Ind is a set of functions from S to R called state indicators (including the constant functions and the set P rop of state propositions, taking values in {0, 1}),

• enabled : S → 2 E are the enabled events in each state

• target : S ×E → S is a partial function describing the state reached from state s on occurrence of an enabled event e ∈ enabled(s)

• delay : S × E → dist(R ≥0 ) is a partial function describing the distribution of the delay of an enabled event e ∈ enabled(s)

• choice : S × 2 E → dist(E) is a partial function describing the distribution to chose among simultaneously scheduled events A configuration of a DESP consists in a state, a value for the current time, and a function E → R ≥0 , called scheduling function, that says for every enabled event when it is scheduled, with value +∞ if an event is not enabled. An evolution step consists in :

• firing the first scheduled event and moving to a new state according to function target, and updating the current time. In case two or more events are scheduled at the same time, select one using the choice distribution.

• update the scheduling function by removing events no longer enabled and sampling a delay for the newly enabled ones

For the sake of saving space in this paper, we omit the formal definition of DESP semantics that can be found in [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] and give an informal description of non Markovian stochastic Petri nets, the high-level language adopted to characterize DESPs in the context of HASL model checking.

DESP in terms of non Markovian stochastic Petri nets.

According to its definition, the characterization of a DESP is a rather unpractical one, requiring an explicit listing of all of its elements (i.e. states, transitions, delay distributions, probability distributions governing concurrent events). However, several high-level formalisms commonly used for representing Markov chain models (e.g. generalized stochastic Petri nets [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF], Stochastic Process Algebras [START_REF] Gorrieri | Unified specification and performance evaluation using stochastic process algebras[END_REF]), can straightforwardly be adapted to represent a quite expressive class of DESPs. In the context of HASL model checking, we consider non Markovian stochastic Petri nets (NMSPNs) [START_REF] Bobbio | Recent developments in non-Markovian stochastic Petri nets[END_REF] as a high level formalism for representing DESPs. Our choice of NMSPNs has been guided in particular by the following two factors: (1) they allow for efficient path generation (due the simplicity of the firing rule which drives their dynamics) and (2) the distributions used to sample firing times of transitions are not restricted to (a small set containing) exponential distributions and can model non-Markovian behaviors.

We quickly recall the basics about NMSPN models pointing out the correspondence with the various parts of a DESP.

A NMSPN model (e.g. Figure 1) is a bi-partite graph consisting of two classes of nodes: places (represented by circles) and transitions (represented by bars) as well as directed edges, called arcs, possibly annotated by a number: its multiplicity. Places may contain tokens (e.g. representing the number of molecules of a given species) while transitions (i.e. representing the events) indicate how tokens "flow" within the net. The state of a NMSPN consists of a marking, i.e. a place-vector indicating the distribution of tokens throughout the places. Given a transition t and a place p, we say that p is an input (resp. output) place of t whenever there is an arc from p to t (resp. from t to p). A transition is enabled whenever all of its input places contain a number of tokens greater than or equal to the multiplicity of the corresponding (input) arcs. An enabled transition may fire consuming tokens (in a number indicated by the multiplicity of the corresponding input arcs) from all of its input places and producing tokens (in a number indicated by the multiplicity of the corresponding output arcs) in all of its output places. Transitions can be either timed (denoted by empty bars, if exponential, or gray bars if non-exponential) or immediate (denoted by black filled-in bars). Generally speaking, transitions are characterized by: (1) a distribution which randomly determines the delay before firing it (corresponding to the DESP delay function); (2) a priority which deterministically selects, among the transitions scheduled the soonest, the one to be fired; (3) a weight, that is used in the random choice between transitions scheduled the soonest with the same highest priority (priority and weight corresponding to the DESP choice function). With the GSPN formalism [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF], the delay of timed transitions is assumed exponentially distributed, whereas with NMSPN, it can be given by any distribution. Thus when a GSPN timedtransition is characterized simply by its weight t ≡ w (w ∈ R ≥0 indicating an Exp(w) distributed delay), a NMSPN timed-transition is characterized by a triple: t ≡ (Dist-t, Dist-p, w), where Dist-t indicates the type of distribution (e.g. Unif), dist-p indicates the parameters of the distribution (e.g [α, β]) and w ∈ R ≥0 is used to probabilistically choose between transitions occurring with equal delay 1 . If the precise description of the chemical process is of no use here, it is worth 1 equal delay occurring with non null probability in case of non-continuous distributions describing the reactions and the net a bit further. One place has been created for each species involved in reactions R 1 and R 2 (i.e. P ro, RNAp, P rox, RBS, RNA), plus an extra place (i.e. RNA) for capturing the intermediate delayed phase of RNA formation. Then reaction R 1 , corresponding to transition init in the net, follows a classical exponential distribution, with in this case rate k t = 0.25. Reaction R 2 is a bit more complex since not all products are released at the same time. After an exponentially distributed delay of rate λ 1 = 1/400, products Pro and RBS are released (transition transc) while RNA needs another Gamma distributed delay (transition termin with parameters shape = 1000 and scale = 0.09 as from experimental data) to be released. In the initial marking M 0 = (1, 2, 0, 0, 0, 0) we assume a molecule of Pro and two molecules of RNAp are available. Thus in state M 0 init is the only reaction enabled, and when it fires it will remove one token from both P ro and RNAp and add a token in each of its output places (i.e. RNAp and P rox), changing the marking into M 1 = (0, 2, 1, 0, 0, 0) whereby the only enabled transition is transc. Upon firing, this transition adds a token back to place P ro which will enable again transition init and lead to a possible alternation of these two transitions.

Pro

Hybrid Automata Stochastic Language

HASL is a formalism designed to analyze properties of a DESP D. It is quite expressive since it permits both to estimate the occurrence probability for a given set of paths and to evaluate complex measures on the system. To do so, a HASL specification is a pair (A, Z) where A is Linear Hybrid Automaton (i.e. a restriction of hybrid automata [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF]) and Z is an expression, involving data variables of A, that states the measure under study. The goal of HASL statistical model checking is to estimate the value of Z by synchronization of the process D with the automaton A. More precisely, a timed execution of D is sampled following the probability distributions. By construction of the LHA, at most one execution of A can follow this execution of D, and this goes on until A reaches some final state or the synchronization fails, which occurs when a transition of D cannot be followed by A. The role of the LHA is twofold. First it permits to select paths: if a path of the corresponding DESP is relevant, the automaton reaches a final location, and if the automaton decides that a path is irrelevant (failed synchronization), it is discarded and not taken into account to compute expression Z. Second the evolution of data variables maintained by the automaton are not only used to select paths but also to get precise information about the system (such as time spent in a particular location, number of molecules of a given type, cumulated waiting time of products in a queue,... ), information that will be used in the final expression Z.

Synchronized Linear Hybrid Automata

In order to evaluate sophisticated measures on the DESP, we need an expressive model for selecting relevant paths while computing interesting information along the paths. To do so we need more than traditional timed automata [START_REF] Alur | A theory of timed automata[END_REF]. First in our model we want to consider not just clock variables but rather a set X = (x 1 , ...x n ) of real valued data variables, whose evolution rate can change depending on the location of the automaton but also on the state of the DESP. Then we also want to extend constraints (the conditions under which it is possible to fire a transition) and updates (the effect on variables of firing a transition). We thus introduce a set Const of boolean combinations of inequalities of the form 1≤i≤n α i x i + c ≺ 0 where α i , c ∈ Ind are indicators and ≺∈ {=, <, >, ≤, ≥}, and a set Up where an update in Up is an n-tuple of functions u 1 , ..., u n in which each u k is of the form x k = 1≤i≤n α i x i + c where the α i and c are indicators.

We then obtain a subclass of hybrid automata, enriched with indicators, that we call linear hybrid automata (LHA), which is formally defined as follows:

Definition 2.2. A (synchronized) LHA is a tuple A = E, L, Λ, Init, Final , X, flow, → composed of:
• E, a finite alphabet of events;

• L, a finite set of locations;

• Λ : L → P rop, a location labeling function;

• Init, a subset of L called the initial locations;

• Final , a subset of L called the final locations;

• X = (x 1 , ...x n ) a n-tuple of data variables;
• flow : L → Ind n , a function that associates to each location a n-tuple of indicators representing the rate of evolution of each data variable.

• →⊆ L × 2 E ∪ { } × Const × Up × L, the transition relation
An edge (l, E , γ, U, l ) ∈→ (also denoted l E ,γ,U ----→ l ), consists of: two locations l and l (source and destination), either a set E of labels telling which events of the DESP can trigger this edge of the LHA, or the label denoting an autonomous (i.e. not synchronized) edge, a constraint γ ∈ Const indicating for which values of the data variables this edge can be fired, and a set U of updates saying how the variables' values are affected on traversing this edge. Updates and constraints are here much more expressive than those authorized in timed automata as they allow for linear combinations of data variables using indicators as coefficients. These LHA are called synchronized since they are meant to synchronize with a DESP through events in E and propositions in P rop. In the following, since all our LHA wil be synchronized, we wil just call them LHA for short. In order to ensure the uniqueness of the path of A corresponding to a run of the DESP, a set of further constraints are imposed on the automaton and can be found in [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF]. For example, locations can be annotated with propositions, ensuring that a location can be entered only when the proposition is true (in the state of the DESP), which in particular can ensure that only one initial location is possible for each initial state of the DESP. Another condition is to restrict autonomous transitions to a set lConst of left closed constraints as they have to be taken "as soon as enabled". Example 2.2. : Figure 2 depicts two variants of a simple, two locations (l 0 initial, l 1 final), LHA for selecting paths of the NMSPN toy model of Figure 1. Such LHA employs two data-variables: t ∈ R ≥0 , registering the simulationtime (i.e. flow ṫ = 1 in every location), and n 1 ∈ N, counting the occurrences of transition transcr (i.e. flow ṅ1 = 0 in every location). The automaton synchronizes with events occurrences (i.e. firing of NMSPN transitions) and distinguishes between the occurrences of the transc event (captured by the top self-loop on l 0 through which n 1 is incremented) and the occurrences of any other event (captured by the bottom self-loop on l 0 ). The autonomous edge l 0 → l 1 (leading to accepting location l 1 ) characterizes the selected paths. Thus the LHA in Figure 2(a) selects time-bounded paths corresponding to simulation time t = T , while the LHA in Figure 2(b) selects event-bounded paths which contains n 1 = N occurrences of event transc.

HASL expressions

The second component of an HASL specification is an expression, denoted Z that is either a particular operator (probability, CDF, PDF, as described later) or the expectation (E) of an arithmetic expression whose basic ingredients are data variables and that uses path operators like last(y) (i.e. the last value of y along a synchronizing path, when reaching a final location), min(y) (max(y)) the minimum (maximum), value of y along a synchronizing path), int(y) (i.e. the integral over time along a path) and avg(y) (the average value of y along a path). • φ a1 ≡ (A a , E[last(n 1 )]): expected number of transc events within time T .

• φ b1 ≡ (A b , E[last(t)]): expected delay for N occurrences of transc.

Measures of probability in HASL terms. Given a set of executions Ex of the DESP, for example the set of executions satisfying a given property, one may want to compute the probability of this set of executions. Measures of probability can be obtained with HASL by building an automaton that only accepts executions in Ex and using a dedicated HASL expression Z ≡ P , the evaluation of which results in the probability that a path of the model is accepted by the LHA (this is simply given by the ratio of accepted paths over the whole set sampled from the model).

Figure 3 shows another variant of the previous LHA, which we may refer to as A c . In this case the selected paths are distinguished with respect to the occurrences of the transc event they contain: for those which contain exactly N occurrences of transc the variable OK is set to 1, for the remaining ones OK is set to 0. Thus an example of HASL encoded measure of probability is:

φ c1 ≡ (A c , P ) : probability of N occurrences of transc within time T
In recent updates, the statistical model checker for HASL (introduced in next section) has been enriched with operators for automatically measuring the Probability/Cumulative Distribution Function (PDF/CDF) of the value that a variable x takes at the end of a synchronizing path. This is essentially achieved by specification of a partition of the domain of the variable x obtained through the following syntax: Z = P DF (last(x), lower, step, stop) respectively Z = CDF (last(x), lower, step, stop), where lower denotes the lower bound of the domain, stop the lower bound of last interval, and step the dimension of each sub-interval. Thus, for example, the expression Z = P DF (last(n 1 ), 0, 1, 10) splits the domain of n 1 ∈ N, into {0}, {1} . . . {9}, {10, . . . ∞}), through which the PDF of the number of occurred transc event within T is sampled at values 0 to 10.

The COSMOS tool: HASL statistical model checking scheme

Assessment of HASL specifications against a DESP model is performed by means of the COSMOS [START_REF]Cosmos home[END_REF][START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] statistical model checker. COSMOS employs a classical confidence interval method [START_REF] Ross | Simulation[END_REF] to obtain an estimate Z of the target expression Z. The running scheme of COSMOS is as follows: given as inputs a NMSPN model D, an HASL specification (A, Z), a confidence-level and a semi-interval width δ, COSMOS applies stochastic simulation to iteratively generate trajectories of the product-process (D × A) in an amount sufficient so that the resulting estimate Z is guaranteed to meet the required accuracy (i.e with probability (1 -) the actual value of Z will fall in the estimated interval whose semi-width is bounded by δ, i.e. [ Zδ, Z + δ]).

Following the characteristics of the above presented models, the synchronization of a NMSPN and an LHA admits two possible behaviors. Either it ends up in some absorbing configuration, corresponding to either a final state of the automaton or a failure of synchronization and leading to a finite synchronizing path, or the synchronization goes over all states of the path without ever reaching an absorbing configuration. In order to perform statistical verification, all the sampled paths need to be finite. In order to ensure such a property, we assume that, with probability 1, the synchronizing path generated reaches a final state. This semantical assumption can be ensured by structural properties of the automaton and/or the NMSPN. For instance, the time bounded Until of CSL guarantees this property. As a second example, the time unbounded Until of CSL also guarantees this property when applied on finite CTMCs where all terminal strongly connected components of the chain include a state that fulfills the target sub-formula of the Until operator. This (still open) issue is also addressed in [START_REF] He | A bounded statistical approach for model checking of unbounded until properties[END_REF][START_REF] Sen | VESTA: A statistical model-checker and analyzer for probabilistic systems[END_REF]. Due to this assumption, the random path variables are well defined and the expression associated to the specification may be evaluated with expectations defined w.r.t. the distribution of a random path conditioned by acceptance of the path. In other words, the LHA both calculates the relevant measures during the execution and selects the relevant executions for computing the expectations.

For a reader interested in statistical model checking, we mention several other tools dedicated to, or including a bit of, statistical model checking. This list is not exhaustive. The probabilistic model checker PRISM [START_REF] Kwiatkowska | Prism 4.0: Verification of probabilistic real-time systems[END_REF], originally only using numerical techniques now includes the possibility to check probabilities using statistical techniques. The well established timed model checker UppAal now has a statistical engine [START_REF] David | Time for statistical model checking of real-time systems[END_REF] in which timed systems are given a stochastic semantics and simulation is used to verify stochastic properties. The model checker Ymer [START_REF] Younes | Ymer: A statistical model checker[END_REF] is fully dedicated to statistical model checking and can handle GSMPs as well as Markovian models, checking CSL properties. We can also mention VESTA [START_REF] Sen | VESTA: A statistical model-checker and analyzer for probabilistic systems[END_REF], a statistical model checker that also computes expected values of several performance evaluation expressions, and PLASMA [START_REF] Jégourel | A platform for high performance statistical model checking -plasma[END_REF], also a statistical model checking platform that verifies (a variant of) bounded LTL properties.

Application of HASL to the analysis of stochastic oscillations

Oscillatory trends are a fundamental aspect of the dynamics of many biological mechanisms, therefore the ability to detect/measure them is crucial. In a stochastic modeling context, the aim is to be able to measure relevant aspects of oscillations such as the average amplitude/period of oscillations or, even better, the probability density function of oscillations' amplitude and period. In this section we discuss the application of the HASL method to measuring these aspects. As a first step, we look at the formal definition of what an oscillation is and provide two distinct characterizations upon which we base the construction of HASL specifications for assessing different aspects (i.e. period duration and amplitude) of oscillatory traces. We then demonstrate the application of such HASL specifications on examples of stochastic oscillators in Section 3.2.

Measuring oscillations with HASL

Among various performance evaluation measures, the HASL formalism can be used for characterizing and effectively assessing oscillation related measures of DESP models. To this aim we consider two orthogonal approaches, i.e. two different automata-based characterizations of oscillatory traces. In the first one we introduce an LHA, A peaks , for assessing the peaks of so-called noisy alternating traces: i.e. traces exhibiting a (possibly infinite) sequence of local maxima interleaved with local-minima. The A peaks automaton allows the user to assess different measures related to detected peaks such as, for example, the average value or the PDF/CDF of the peaks of a trace. In the second approach (based on [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]), we introduce an LHA, denoted A per , for assessing the duration of so-called noisy period for noisy-periodic traces. i.e. traces that perpetually switch between states of a low to states of a high partition of the population domain of the observed species. The A per automaton allows the user to assess different measures related to the periodicity of an oscillator such as, for example, the average duration as well as the fluctuation of the period of a stochastic oscillator. As we will see, both automata depend on certain configuration parameters and the chosen configuration establishes an observational perspective over the oscillating traces. For example with A per the considered high/low partition of the observed species state space affects the measured period. In this respect we suggest a two-step approach should be followed: given an oscillating model, first the A peaks LHA is applied, resulting in estimates of the peaks (hence the amplitude) of oscillations. Peaks measurements can then be used to identify a sensible partition of the domain of the oscillating species, hence to properly set up the parameters of the LHA for measuring noisy periodic traces. Figure 13 illustrates for example how, using automaton A peaks , it is possible to set a H threshold capturing 95% of the local maxima. An automaton using these thresholds will then be used to assess the duration of oscillation periods. We give details of both methods in the following.

Notation and vocabulary. Referring to the characterization of the traces of a DESP model, we introduce few basic notations and concepts we will refer throughout the remainder of this section. We denote σ an (infinitely long) simulation trace of an n-dimensional DESP model whose states' form is s = (s 1 , . . . s n ) ∈ N n , with s i being the value along the ith dimension (e.g the number of molecules of species i). We denote σ i (t) the i th projection of σ(t), the state in which σ is at time t. Following the characterization given in [START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in Markovian population models[END_REF], we say that σ i is either: convergent (i.e. tending to a finite value), divergent (i.e. tending to infinity) or oscillating (i.e. the lack of the previous two). Furthermore, σ i is periodic with period δ > 0 iff ∀t, σ i (t) = σ i (t + δ). Whenever there exists a positive minimal value with this property, it is called basic period. Thus σ is periodic oscillatory along the i-th dimension iff σ i is both oscillating and periodic. We point out that, for stochastic models, "exact" periodicity (as given above) is too restrictive hence in the following (Section 3.1.2) we introduce the less restrictive definition of noisy periodicity [START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in Markovian population models[END_REF].

Measuring noisy alternating traces

We consider the problem of detecting the local maxima (minima) on traces of a given DESP model. 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. In the most general case every single point corresponding to a change of trend is accounted as a maximum/minimum, even if distanced of only one unit from the previous one. In this case, however, the detected maxima (minima) will include the minimal peaks corresponding to stochastic noise, hence actually biasing local maxima (minima) related measures with the effect of stochastic noise. As a consequence, we propose a parametric characterization (Def. 3.1) of local maxima (minima) based on a so-called noise parameter which establishes an observational perspective: a point corresponding to a change of trend is detected as local maximum (minimum) only if its distance from its predecessor is at least equal to the value of the noise parameter. We then introduce the automaton A peaks which encodes such a characterization and provides us with an effective machinery to detect local maxima/minima over traces of a DESP model. Definition 3.1. Let σ be a trace of an n-dimensional (n ≥ 1), possibly infinitestate, DESP, and let δ ∈ N ≥0 represent a noise level. Trace σ i is said noisy periodically alternating (or noisy alternating for short), with respect to noise level δ, iff it perpetually alternates from a local minimum to a local maximum and the minimal distance between consecutive minima and maxima is no smaller than δ.

The detection of local maxima/minima along an alternating trace is illustrated in Figure 4. The A peaks automaton: We introduce an LHA (Figure 5), denoted A peaks , designed for detecting local maxima/minima along noisy alternating traces of

start ṫ : 1 noisyDec ṫ : 1 noisyInc ṫ : 1 Min ṫ : 1 ] ,( A < x ) ^(n M < N 1 ) , { n M + + , L m a x [x ] + + , S m a x + = x , x := A } ] , ( A < x ) , { x : = A } E + A , > , ; E=A,>,; Max ṫ : 1 E+A,>,{x := A} E=A,>,; ] , ( A > x ) , { x : = A } E,>,; E,>,; ], (A > x + ), { x := A } ] ,( A < x ) ,{ 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 a given observed species (in this case A). The automaton A peaks is parametric with respect to a parameter δ. This parameter represents the level of noise tolerated in the system: after a local maximum, if the population decreases of less than δ, we keep looking for a higher local maximum. If the decrease is greater than δ then the maximum has been found and the automaton now looks for the next local minimum. Given the observed quantity A, the definition of A peaks requires partitioning the set of events in three subsets:

],(A>x+ )^(nm = N 1), {Lmin[x]++, Smin+= x, nm ++} ],(A<x )^(nM = N 1), {Lmax[x]++, Smax+= x, nM ++} ] ,( A < x + ) ^(n m < N 1 ) , { n m + + , L m i n [x ] + + , S m i n + = x , x := A }
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, registering the minima/maxima along the path. 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 5). During the first initT time units, the simulated trace simply grows without any peak detection2 . The actual analysis begins in location start from which we move to either Max or Min depending whether we initially observe an increase (i.e. N n array of frequency of heights of detected maxima (minima)

Table 1: The data variables of automaton A peaks of Figure 5 for locating the peaks of a noisy oscillatory traces 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 maxima/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 minimum Lmin[x] (maximum Lmax[x]) 3 before storing the new value of A in x and finally increasing 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 a reaction of E +A (E -A ), in which case we also store the new increased (decreased) value of A, hence a potential next local maximum (minimum) in x, or one of E =A . On the other hand, on occurrence of a "decreasing" ("increasing") reaction E -A (E +A ) we

3 with a slight abuse of notation we refer to Lmin[] and Lmax[] as arrays whereas in reality within COSMOS/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 with a negligible 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.

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 exceed x ( x exceed A), or eventually enter Min (Max) as soon as the observed decrease (increase) goes beyond the chosen δ (see above). For the automaton A peaks depicted in Figure 5, 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.

Measures associated with automaton A peaks . We introduce the following four target measures, given in terms of HASL expressions, associated with A peaks . All of them are computed on paths that end as soon as at least N minima or N maxima have been detected, so the measures only concerns transient behaviors up to seeing N extrema of the same type. Observe that a timebounded version of A peaks can be straightforwardly obtained in order to ensure termination (see discussion in Section 3.1.3).

• Z max ≡ E[last(Smax)/n M ]: corresponding to the expected value of the average height along a path of the maximal peaks.

• Z min ≡ E[last(Smin)/n m ]: same as Z max but for minima.

• Z P DF max ≡ P DF (last(Smax)/n M , s, l, h): corresponding to the PDF of the average height (along a path) of the maximal peaks

• Z P DF min ≡ P DF (last(Smin)/n m , s, l, h): same as Z P DF max but for minima

• Z Lmax ≡ E(last(Lmax[k]
)/n M ): expected proportion of maximal peaks of height k along a path. This measure is used to compute the PDF of the height (along a path) of the maximal peaks

• Z Lmin ≡ E(last(Lmin[k]
)/n m ): same as Z Lmax but for minima.

Measuring noisy periodic traces

Here we consider the problem of studying the periodicity of traces of a DESP model. As mentioned earlier, exact periodicity never occurs in practice in non degenerative stochastic systems. Following the approach of [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in Markov population models[END_REF], we introduce a less restrictive characterization of periodicity, namely the so-called noisy periodicity (see Definition 3.2) which is based on two thresholds parameters, L and H. As a consequence, the perceived noisy periods depends on the chosen L and H thresholds. Definition 3.2. Let σ be a trace of an n-dimensional (n ≥ 1), possibly infinitestate, DESP, and let L, H ∈ N (0 ≤ L < H) be a lower, respectively an upper bound, inducing intervals low = (-∞, L], mid = (L, H) and high = [H, +∞), trace σ i is said noisy periodic with minimal amplitude H -L iff it perpetually switches from low to high (possibly passing through mid) and returning to low. Remark 3.1. A trace that is noisy periodic is also necessarily noisy periodically alternating (i.e. noisy-periodicity is a sufficient condition for noisy-periodicalternance), however the opposite is not the case: a periodically alternating trace may well diverge.

According to Definition 3.2 a noisy-periodic trace corresponds to the following ω-regular expression4 :

e np = (low(low + mid) * high(high + mid) * ) ω (1)
Notice that such a characterization induces a bound on the minimal-amplitude of oscillations (corresponding to the width of the chosen partition, i.e. the distance H-L), but entails no specific constraints on the maximal amplitude nor on the period duration of oscillations (i.e. traces belonging to this class will include also those with non-constant period duration, i.e. those with noisy-period). It is also worth mentioning that, since we use an approach based on sampling, we will not be able to study sustained (i.e. perpetual) oscillations and we thus develop different ways to study oscillating characteristics on bounded paths. Figure 6 shows an example of a noisy periodic trace of a DESP. A period corresponds to the behavior of the system between two consecutive crossing points corresponding to entering in one extreme of the partition, interleaved by at least one crossing point into the opposite extreme. For example in Figure 6 the first period p 1 corresponds to the duration between the first mid-to-low crossing and the successive mid-to-low crossing interleaved by a mid-to-high crossing. Notice that the first complete period(p 1 ) might be preceded by a spurious period (i.e. p 0 ) which should be discarded as there's no guarantee that t = 0 corresponds with an actual mid-to-low crossing, hence with the actual beginning of p 0 .

The A per automaton: We introduce an LHA (Figure 7), denoted A per , designed for detecting and measuring the periods of a given observed species (in this case A). More specifically, A per is designed for assessing two characteristics of an oscillatory trace: the mean period duration (denoted tp ) and the period fluctuation (denoted f tp ) over the first N periods detected along a trace. Period fluctuation represents how much (on average) a single period duration deviates from the mean duration computed over the N periods observed along a trace. From the point of view of analysis, the period fluctuation is a useful indication of the irregularity of the observed oscillator. Automaton A per 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). As with A peaks , the simulated trajectory initially unfolds for a given duration (initT ), letting A per within l 0 without doing anything. After initT time units, A per enters low5 (possibly through l 0 if at that time the population is above level L) where the actual oscillation analysis begins. From low the automaton follows the profile of A hence moving to mid as soon as L < A < H holds, and then to high as soon as A ≥ H. With A per , a period starting point is associated with the first mid-to-low crossing that follows a mid-to-high crossing 6 . Hence the first detected period (crossing from low to high and back to low) is discarded as its duration may be spurious. A per uses six variables (Table 2): t is a timer that keeps track of simulation time; n counts the number of detected noisy-periods while top is a boolean flag used for distinguishing the mid to low crossing points that correspond to the closure of one period (i.e. when a name domain

update definition t R ≥0 reset n N increment top bool complement tp R ≥0 reset tp R ≥0 f ( tp, tp, n) = 1 n+1 ( tpn • n + tp) ft p R ≥0 g(ft p , tp, tp, n) = 1 n [(n -1)ft p + (tp -tp)(tp -f ( tp, tp, n + 1))]
Table 2: The data variables of automata Aper of Figure 7 for measures of noisy-periodicity traversal from mid-to-low has been preceded by a traversal from high-to-mid) from those which do not. Notice that, in order to ignore the first potentially spurious period (p 0 in Figure 6), n is initially set to -1 on entering low from the initial transient filter, and the simulation time t is then reset on detection of the closure of the first spurious detected (i.e. on entering low from mid when n = -1). Furthermore t p stores the duration of the last detected period, while t p maintains the mean duration of all (so far) detected periods and f tp stores the fluctuation of the period duration for all (already) detected periods. Notice that the fluctuation (i.e. f tp ) is computed on the fly (see Table 2) by 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. Finally the analysis of simulated trajectories stops (by entering the accepting location end) as soon as the N th period has been detected.

Measures associated with automaton A per . We introduce the following four target measures, given in terms of HASL expressions, associated with A per :

• Z 1 ≡ E[last( tp )]: corresponding to the mean value of the period duration for the first N detected periods.

• Z 2a ≡ 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 2 b ≡ CDF ( tp ), s, l, h): as Z 2a but corresponds to the CDF of the average period duration over the first N detected periods.

• Z 3 ≡ E[last(f tp )]: corresponding to the fluctuation of the period duration.

Expressions Z 1 , Z 2a and Z 2 b represent measures related to the duration of the oscillation periods (i.e. the mean value, and the PDF, respectively the CDF, of the period duration). On the other hand Z 3 is designed to asses the fluctuation of the period duration, 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. Thus an oscillator whose trajectories are likely to exhibit periods of, more or less, the same duration (e.g. Figure 10(a)) will result in smaller fluctuation f tp than an oscillator whose trajectories exhibits periods of more variable duration (e.g. Figure 10(b)). We give an example of analysis of the regularity of the oscillation periods on a model of the repressilator (see Section 3.2).

Remark. As we argued before, the measurements of the period of oscillations obtained through automaton A per are affected by the considered L, H partition, i.e. the chosen value of the L, H parameters of A per . We provide an example illustrating this aspect in Section 3.2.2 where we analyze the repressilator model. The problem is that, without any prior knowledge about the observed oscillator, and in particular about the height of the maximal/minimal peaks of the corresponding traces, we risk to set the L, H parameters to values which are not representative of the actual observed oscillator. For example if we set H to a value which is well above the average height of the maximal peaks then we are going to miss most of the peaks. As a consequence we envisage that, when assessing oscillators, a protocol involving the A peaks and A per automata should be followed. More specifically: first apply A peaks to get a measurement of the distribution of maximal/minimal peaks' values, then use this information to set reasonable values for the L, H parameters of A per . We provide details about the application of such a protocol through an example presented in Section 3.2.1.

Discussion

With respect to the HASL based oscillation analysis approach introduced above, there are a number of delicate aspects which are worth discussing.

Measuring vs verifying. First it is worth stressing that the HASL approach is not a "pure" verification approach, in the sense of classical model checking. A "pure" verification approach would allow to decide whether a given stochastic model oscillates or not. In stochastic modeling deciding about a qualitative property (such as the model oscillates) boils down to comparing a quantitative measure (of probability) against a given threshold. For oscillations this has been shown to be feasible in some cases, for example: whether a finite-state CTMC model oscillates sustainably can be verified [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in Markov population models[END_REF] but requires the computation of the steady-state distribution of (a timed-automaton expanded version of) the considered CTMC 7 (thus is affected by the state-space explosion problem). HASL, as any other statistical model checking approach, is inherently limited to reasoning about finite observations, hence is not suited for verification of infinite behaviors, such as sustained oscillations. Nonetheless, it can be used to extract significant (finite horizon) knowledge which can be used as an indicator of a model's infinite behavior. Thus, if a sustained oscillator is analyzed through HASL, its oscillations will be captured by the A per and A peaks automata whose parameters will be appropriately set so as to discard a possible non-oscillatory initial transient. On the other hand, if a non-oscillating model is studied using the A per and A peaks automata, the number of detected periods will either be null (if the model does not include any period in its initial transient), or it will reach an asymptotic value as we enlarge the observation windows (if the model include some initial periods before either converging or diverging).

On automata's parameters and termination. With respect to termination, observe that A per and A peaks have been, for the sake of simplicity, presented (Figure 7 and Figure 5) as event-bounded measures, thus measuring terminates once a finite number N of oscillation-related events have been observed (i.e. N noisy-periods for A per , N local maxima for A peaks ). This may lead to non-termination (or to a very long runtime) when applied to certain models. To avoid non-termination issues one can simply consider the time-bounded version of both automata, denoted A t≤T per and A t≤T peaks , obtained by adding a time constraint t ≤ T (T being a parameter of the automaton) to each arc of the automaton (hence enforcing the synchronization to terminate, at latest, after T time units).

Another important aspect is the setting of the automata's parameters (L and H for A per , δ for A peaks ) as these establish an observational perspective, hence affect the outcome of the measure. For A per , the further H (resp. L) is chosen above (resp. below) the actual average value of maximal (resp. minimal) peaks of the considered oscillator, the less likely it becomes for a trace to enter the [H, ∞) (resp. [0, L)) region. A simple way to chose "reasonable" L and H values is to plot a few single traces generated for the considered model and to fix the values of L and H close to the peaks as perceived on the traces plot. A more accurate manner to proceed is to first employ A peaks for measuring the PDF of the heights of peaks and then exploit such measure to set A per 's L and H values. Since the outcome of an A peaks measure depends on the chosen δ parameter (i.e. the height of noisy spikes to be filtered out from the detection of local maxima/minima), the measurement of the PDF of peaks should be repeated, starting form δ = 1 and for increasing values of δ. The comparison of results obtained for different small values of δ, (e.g. δ = 1 to δ = 5) will provide information about the localization of noisy spikes, thus allowing for figuring out what is the most likely height of the "actual" peaks of oscillation. In Section 3.2 we illustrate the combined application of A peaks and A per measurements on a concrete example. Initial transient. Another relevant aspect is related to the fact that what we aim at is measuring relevant characteristics of oscillations (e.g. the average duration of the oscillation period) but only once the system has reached its long term behavior, hence disregarding the (spurious) effect that the initial transient might have on the measured outcome. The elimination of the effect of the initial transient from measures obtained via simulation of a stochastic model is a well established subject in the literature [START_REF] Eickhoff | Detecting the duration of initial transient in steady state simulation of arbitrary performance measures[END_REF]. In this respect it would be desirable to equip the COSMOS tool with facilities for automatically detecting the duration of the initial transient and for the elimination of its effects from the target measure. In this paper we only provide a partial (not automatic) scheme to tackle that issue. In particular we equip each LHA dedicated to assessing oscillations with an initial transient filter, i.e. an elementary part of the automaton whose only goal is to let time pass (while simulating a trace) before actually starting the measuring process. Thus the user should first run a number of preliminary experiments (i.e. reiterated with increasing values of the initT parameter) aimed at establishing at which point the estimate of the target measure starts stabilizing . For example our experiments indicated that in the case of the repressilator model of Section 3.2 the initial transient has no effect on the measures of the period as such measures turn out to be independent of initT . A similar question arises with the appropriate length of paths in order to have a measure that is close to what could be computed on infinite paths. Figure 17 shows how, by increasing the length of observed path, we see a stabilization in the measure that gives us a good idea of the long term behavior.

Case study: measuring oscillations of a three-genes repressilator model

To demonstrate the above discussed approach, we have considered different examples of oscillators. As a first step towards the empirical validation of HASL based measurement of oscillations, we have considered a toy model, i.e. a probabilistic square wave oscillator. The goal in that respect was to have a simple model for which an analytical expression, for both the mean value and the fluctuation of the cycle duration, existed. Hence we have compared measurements obtained through A per with the corresponding exact values computed analytically for the square wave oscillator: the estimates obtained with HASL turned out to precisely match the exact values (for details see Appendix A.1). In the remainder of this section, on the other hand, we report on the application of HASL based oscillation analysis to a realistic model of genetic oscillator, known as the repressilator [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF].

A three-genes repressilator

Repressilator model. We consider a model of the three-genes repressilator, a simple genetic network with three elements inhibiting the production of each other in a cyclic mono-directional manner. The repressilator is a real genetic circuit, consisting of the natural repressor proteins and the genes coding for them. The cyclic topology gives a certain robustness to the circuit (Fig. 8) that exhibits oscillations for a broad range of the parameters' values.

21

Here we consider a version of the repressilator model taken from [START_REF] Ribeiro | A general modeling strategy for gene regulatory networks with stochastic dynamics[END_REF] which essentially is a simplified version of the more general model of Elowtitz [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. In our model, gene expression is considered as a single step process occurring when the RNA-polymerase (RNAp) reacts with the promoter (P ro i ). The result of this reaction is the appearance of the newly synthesized protein molecule (P i ) and the release of the reactants, i.e. promoter and RNA-polymerase. The i-th promoter (where i ∈ {1, 2, 3}) yields the corresponding protein (P i ), and this is modeled by reaction [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF].

P ro i + RNAp k -→ RNAp + P ro i + P i (2)
Notice that in our model RNAp is never consumed hence its only role is to control the speed of the expression of the three genes (higher RNAp population induces faster gene expression). The repression reactions are modeled as binding/unbinding of the repressor protein to/from the corresponding promoter providing the repressilator's cyclic topology. These reactions are shown in [START_REF] Alur | A theory of timed automata[END_REF]. Protein (P j ) binds to the promoter (P ro i ), forming the complex (P j P ro i ) which prevent the further expression reaction for promoter P ro i . The unbinding of the protein makes the promoter available for the expression again. In reactions (3), i ∈ {1, 2, 3} and j = i -1, except for i = 1, in which case j = 3.

P j + P ro i kr -→ P j P ro i P j P ro i ku -→ P j + P ro i (3)
Finally, proteins undergo degradation (reactions (4)), which is modeled as a single step uni-molecular reaction. Protein degradation takes place either when proteins are free (with kinetic rate d), or when they are bound to promoters (with kinetic rate d ). The reactions (4) take the same values for i and j as in reactions [START_REF] Alur | A theory of timed automata[END_REF].

P j d -→ ∅ P j P ro i d -→ P ro i (4)
Within the stochastic interpretation, the reactions (2), ( 3) and ( 4) are associated to (negative) exponential distributions, hence the repressilator can be simulated following Gillespie's [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] approach. In the remainder we will mostly consider the following default configuration of the model parameters: k = 0.01, k r = 1, k u = 0.08, d = d = 0.01, with RNAp = 30 being the initial population of RNA-polymerases and also assuming that the system contains only one copy of each gene, i.e. every promoter is present in a single copy. However we will also analyze the model's dynamics under perturbation of the following parameters: the RNAp initial population, and the degradation rate (d ) of proteins when attached to the promoter regions. The complete set of 15 reactions for the three-genes repressilator model is given by equations (5), while its encoding in GSPN form, which is what we used for our experiments with the COSMOS tool, is illustrated in Figure 9. Figure 10 shows the dynamics of the three expressed proteins along a single simulation trace of the repressilator model referred to different initial populations of the RNAp (i.e RNAp = 30 and RNAp = 200). The three proteins' populations oscillate in a co-ordinated manner, however the curves are rather noisy (with respect to both amplitude and period duration). Observe that increasing the RNAp population, hence the speed of protein expression, induces a much more pronounced irregularity of the period duration (Figure 10(b)).

R3 1 :P1P ro2 ku --→ P1 + P ro2 R4 1 :P1P ro2 d -→ P ro2 R3 2 :P2P ro3 ku --→ P2 + P ro3 R4 2 :P2P ro3 d -→ P ro3 R3 3 :P3P ro1 ku --→ P3 + P ro1 R4 3 :P3P ro1 d -→ P ro1 R5 1 :P1 d -→ ∅ R5 2 :P2 d -→ ∅ R5 3 :P3 d -→ ∅ (5) 30 RNAp Pro1 Pro2 Pro3 P1 P2 P3 P1Pro2 P2Pro3 P3Pro1 R1 1 R2 1 R1 2 R1 3 R2 2 R2 3 R3 1 R3 2 R3 3 R4 3 R4 2 R4 1 R5 1 R5 2 R5 3
Strategies for the analysis of the repressilator model. As in most realistic stochastic oscillators, oscillation traces in the repressilator are very noisy (see Figure 10), and dependent on several parameters of the model (e.g. initial populations, decay rates, etc.). In Section 3.1.2 we have introduced the A per automaton (Figure 7) as a means for measuring noisy periodic traces. However, the resulting measures are affected by the chosen L and H threshold parameters (i.e. intuitively: the larger the distance H -L the longer the period duration). A simple manner to find out reasonable values for the L and H parameters consists in generating a simulation trace and manually establishing how to partition the Figure 10: A single simulation trace of the three-genes repressilator model with initial RNAp =30 (left) and RNAp =200 (right): observe that with higher initial RNAp the trace consists of higher peaks grouped into rather irregularly distanced batches, hence, in general, higher RNAp population seems to result in increased irregular oscillations.

domain of the oscillating species so that the high (resp. low) interval contains the maximum (resp. minimum) peaks that we consider as those characterizing a period. Despite its simplicity, this approach has important drawbacks: first it is not formal (establishing the L and H values by observation of a single trace is a poorly accurate approach, with no control on the introduced approximation); second it is not automatized.

In the following we propose a methodology for establishing reasonable L and H values based on the measurement of the PDF of local maxima (minima) of the repressilator.

Detecting the most-likely local-maxima (minima) of the repressilator. Figure 11 and Figure 12 depict the PDF of the local maxima, respectively minima, of the repressilator measured through the A peaks automaton (Figure 5). Such measures have been done for different values of the proteins' (i.e. P 1 , P 2 , P 3 ) decay rate: d = 0.01 (the "standard" decay) and d = 0.02 (a two times faster decay). To assess the effect that the chosen level of noise (i.e. the δ parameter of A peaks has on the the outcome, we repeated the experiment for different values of A peaks 's noise parameter, hence both Figure 11 and Figure 12 contain several plots corresponding with different values of noise.

For the PDF of local maxima given in Figure 11, we can see that all plots exhibit a truncated bell-shape profile with truncation, as expected 8 , corresponding to the chosen level of noise. Furthermore, for small values of noise the probability mass exhibits a truncation peak centred on the truncation point. Such a peak reflects the specific nature of the oscillations of the repressilator, and more specifically it provides insights about where noisy small spikes are statistically located. By looking at single traces of the repressilator (e.g. Figure 10(a)), we observe that noisy small spikes appear to be more frequent around the minimal and maximal peaks of the oscillations than along the steep parts, thus when noise is set to small values such small spikes are also detected by A peaks and this yields the spurious peak on the PDF.

Observe further that for sufficiently large values of noise (i.e. noise ≥ 5) the PDF plots are equally centered on a value (≈ 15, for decay rate = 0.02, and ≈ 26, for decay rate = 0.01) which corresponds with the average of the detected local maxima (see Figure 14(a)).

For the PDF of local minima shown in Figure 12, (almost) all plots exhibit a monotonically decreasing profile with a common maximum at zero, which indicates that the most likely minimum of each oscillation period is indeed located at zero (in line with the profile of single traces such as e.g. Figure 10). However, for low values of noise, the probability mass (Figure 12(d) and Figure 12(c)) actually includes a little bump slightly below the region where the local maxima probability mass is concentrated (≈ 15, for decay rate = 0.02, and ≈ 26, for decay rate = 0.01). Again such bump is a reflection of the stochastic character of the repressilator and indicates that little noisy spikes are also more densely located towards the upper peaks of oscillations. Observe that such bumps correctly fades away as we increase the level of noise in our experiments (i.e. we filter out the noisy spikes from the measure).

Exploiting the PDF of maxima/minima to set the L,H thresholds. An issue with the A per automaton is that it requires to set the L,H threshold parameters, whose value has also an effect on the measured periods (see Figure 17 maxima and minima estimated with the A peaks automaton. Intuitively for a given probability value p (i.e. quantile) we chose H (L) so that the probability that local maxima (minima) are above (below) H (L) is at least p. This coincides with determining the quantile of the CDFs of local maxima (minima). Figure 13 shows the reverted CDFs (i.e. (1 -CDF )) of local maxima (straightforwardly derived from the PDFs in Figure 11). By choosing p ∈ (0, 1) and looking at the x-coordinate of the intersection points of each curve with y = p we obtain the desired value for H. For example, by setting H = 13 (x-coordinate of the intersection of yellow curve with p = 0.95) we know that with 95% chance, when protein P 1 decays with rate 0.01, the local maxima of oscillations of minimal amplitude 10 (i.e. noise = 10) will be above H.

Studying the effect of protein degradation on the repressilator oscillations

By application of automata A per and A peaks , we have analyzed the effect that the decay rate (denoted d 1 hereinafter) of protein P 1 has on the repressilator oscillations. Speeding up d 1 corresponds to augmenting the throughput of P 1 degradation, thus intuitively should result in a lower P 1 average population, hence a lower height of maximal peaks of P 1 oscillation, and, as a consequence, 1 CDF (local maxima) X X X X X X 10.5

Figure 13: Setting the H threshold parameter with respect to the 95% quantiles of the CDF of the local maxima of the repressilator. On the x axis, in red, the values above which the local maxima of oscillations fall with 95% chance (e.g. with decay rate = 0.01 the 95% of the local maxima, which are at least noise = 10 units distant from the preceding local minimum, will be greater or equal to 13. Thus, by setting H = 13, we are guaranteed to miss out no more than 5% of all observed such local maxima.) also a shorter period duration (i.e. the delay between two peaks should shorten as the height of peaks decreases). Conversely, slowing down d 1 should have the opposite effect (i.e. a higher average population, hence higher oscillation peaks, hence longer period duration). Figure 14 As for the period duration (Figure 14(b)), we can notice the effect that the observational perspective, i.e. the chosen L, H parameters of A per , has on the target measure, i.e. the average period duration. The average period measured for trajectories with maximal peaks above 8 and 10 (i.e. curves H=8 and H=10) monotonically decreases over d 1 ∈ [0.005, 0.02] thus confirming the intuition that the period of oscillation shrinks as we speed up P 1 's degradation. However, if we increase the threshold and only consider trajectories with peaks above H = 15, the measured period turns out not to be monotonically decreasing but rather displays a minimum at approximately d 1 = 0.015. This, which may look as an anomaly, is actually a sensible consequence of the configuration of the H parameter of A per . What happens in reality to the P 1 signal as we increase d 1 is that the peaks of oscillations are loosing height hence they are less likely to go above a H = 15. On the hand with H = 15 the period is computed only based on these (less likely) high peaks. This artificially increases the computed period.

Assessing the robustness of the repressilator oscillations

As an application of the measurement of the fluctuation of the cycle duration using HASL, we considered the problem of analyzing how robust the oscillations of the repressilator are with respect to certain parameters of the model, in particular: 1) the rate at which proteins are expressed (which in our model is proportional to the population of RNAp) and 2) the rate d at which proteins degrade while attached to the promoter of the gene they repress (i.e. the rate of the reactions P ro j P i d -→ P ro j ). We have already pointed out (Figure 10) that an increase of the protein expression speed (i.e. through RNAp population) induces an increased irregularity of the periodicity of the peaks. Such an irregularity is explained as follows: with faster expression proteins have a tendency to rapidly accumulate (in places P 1 , P 2 , P 3 ), hence increasing the repression of the controlled protein by constantly occupying the corresponding promoter (higher throughput of reactions R 2x ). One way for compensating the effect of faster expression is simply to speed up protein degradation, for example through rate d . Figure 15 shows traces for configurations corresponding to a factor 33 increase of the expression speed (i.e. RNAp = 1000) and for different values of degradation rate d . In absence of compensation (d = 0.01 as in default configuration) a 33x speedup of the expression process results into a degenerated signal consisting of irregularly scattered spikes, superposed to a noisy base (Figure 15(a)). However, as we increase d we progressively re-establish the oscillatory regime of the repressilator (Figure 15(d)).

If the traces in Figure 15 illustrate experimentally the (ir)regularity of oscillations, it is possible to measure this irregularity using HASL, as we will explain later. Figure 16 shows the value of the average cycle duration and the average cycle fluctuation, measured (through automaton A per , expressions Z 1 , resp. Z 3 , and using L = 1, H = 100 as thresholds for detecting periods) for the speeded up (i.e. RNAp =1000) model and for different values of the degrada- 15 where we move from trajectories consisting of irregularly batched peaks (e.g. Figure 15(b)) towards trajectories consisting of more evenly distributed peaks (e.g. Figure 15(d)).

On the convergence of estimates in function of the observation length

We stress that the HASL based measures of oscillations discussed in this paper are, inherently, approximations of the actual quantities (e.g. average cycle duration and fluctuation) of a stochastic oscillator model. This is due to the fact that HASL is a statistical, rather than a numerical, model checking approach, therefore all quantities are calculated with respect to (sampled) trajectories of finite length (i.e. either time-bounded or event-bounded trajectories). In order to assess how fast the measured quantities (cycle durations and peaks of oscillation) converge for the repressilator model we performed a number of re-iterated experiments where we evaluated the same quantity but changed the length of each observation (i.e. the number of observed cycles or peaks in each trajectory). Figure 17 17(a) show that the measured cycle duration increases as we increase H. This result is not surprising as, by increasing H (while keeping L constant), we essentially instruct our cycle detector to exclude from detection those maximal peaks of lower height, hence, unless we change the stochastic characteristics of the oscillator, this necessarily implies that it will take longer for a trajectory to reach a higher H crossing point. Similarly Figure 17(b) reports on re-iterated measures of the min/max peaks obtained by using different values of the noise parameter δ of A peaks . Observe that in Figure 17(b) the values of maximal peaks are represented by full lines (each color corresponding to a different value for the noise parameter δ) while that of corresponding minimal peaks by a dashed line of the same color.

Concerning the minimal peaks, the mean value is between 0 and 1 for noise 10 to 15 whereas it jumps to 6 for noise level 5. This suggests that a noise level of 5 is too small, hence detecting spurious minimal peaks.

Related work

The application of temporal logic reasoning to the analysis of oscillatory trends of stochastic models has been considered before. CSL-based characterizations of oscillations for CTMC models of biochemical reactions have 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 | Characterizing oscillatory and noisy periodic behavior in Markov population models[END_REF][START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in Markovian population models[END_REF]. In particular in [START_REF] Spieler | Characterizing oscillatory and noisy periodic behavior in Markov population models[END_REF] it has been shown that oscillations can be detected with CSL through an automata-based procedure. This, however, requires a manual hardwiring of the detector-automaton within the CTMC model. The procedure has been demonstrated by implementation of examples on the PRISM model-checker. The main drawback of such an approach is indeed due to the required hardwiring of the automaton within the CTMC model, a rather costly practice from the modeling point of view, which provides little flexibility. Recently, measuring of oscillations have been considered with other statistical model checking tools, i.e. UPPAAL-SMC [START_REF]Uppaal-smc home page[END_REF], using MITL [START_REF] David | Runtime verification of biological systems[END_REF] and PLASMA [START_REF] Jégourel | A platform for high performance statistical model checking -plasma[END_REF], using BLTL [START_REF] Finkbeiner | Checking finite traces using alternating automata[END_REF]. However, the approach used for detection of the cycle duration in that case is inherently limited: i.e. the level of nesting of a MITL formula for detecting multiple occurrences of cycles is proportional to the number of cycles one wants to detect. More generally, the approaches in [START_REF] Ballarini | Query-based verification of qualitative trends and oscillations in biochemical systems[END_REF][START_REF] David | Runtime verification of biological systems[END_REF] look at a simplification of the much more generic problem of oscillation detection whose complexity has been nicely illustrated in [START_REF] Spieler | Model checking of oscillatory and noisy periodic behavior in Markovian population models[END_REF]. Another interesting contribution in the field of analysis of oscillations 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 analyzed, however by mathematical approaches rather than temporal logic based ones.

Analysis of Single Gene model with stochastic delayed dynamics

We illustrate the application of HASL model checking to the analysis of a model of gene expression with stochastic delayed dynamics. We (1) present the NMSPN encoding of the considered model; (2) introduce a number of relevant properties/measures of a model, first describing them informally and then providing their encoding in HASL terms; (3) discuss results obtained by evaluation of the presented properties/measures by means of the COSMOS model-checker.

Stochastic models of gene expression with delayed dynamics

Gene expression is the process by which proteins are synthesized from a sequence in the DNA. It consists of two main phases: transcription and translation. Transcription is the copying of a sequence in the DNA strand by an RNA-polymerase (RNAp) into an RNA molecule. This process consists of three main stages: initiation, elongation and termination. Initiation consists of the binding of the RNAp to a promoter (Pro) region, unwinding the DNA and promoter escape. Afterwards, elongation takes place, during which the RNA sequence is formed, following the DNA code. Once the termination sequence is reached, both the RNAp and the RNA are released. In prokaryotes, translation, the process by which proteins are synthesized from the (transcribed) RNA sequence, can start as soon as the Ribosome Binding Site (RBS) region of the RNA is formed. The rate of expression of a gene is usually regulated at the stage of transcription, by activator/repressor molecules that can bind to the operator sites (generally located at the promoter region of the gene) and then promote/inhibit transcription initiation. Evidence suggests that this is a highly stochastic process (see, e.g. [START_REF] Arkin | Stochastic Kinetic Analysis of a Developmental Pathway Bifurcation in Phage-l Escherichia coli[END_REF]) since, usually, the number of molecules involved, e.g. transcription factors and promoter regions, is very small, ranging from one to a few at a given moment [START_REF] Taniguchi | Quantifying E. coli Proteome and Transcriptome with Single-Molecule Sensitivity in Single Cells[END_REF]. Consequently, stochastic modeling approaches were found to be more appropriate than other strategies (e.g. ODE models or Boolean logic [START_REF] Szallasi | System Modeling in Cellular Biology: From Concepts to Nuts and Bolts. A Bradford book[END_REF]).

The first stochastic models of gene expression assumed that the process of gene expression, once initialized, is instantaneous [START_REF] Arkin | Stochastic Kinetic Analysis of a Developmental Pathway Bifurcation in Phage-l Escherichia coli[END_REF]. Namely, each step was modeled as a uni-or bi-molecular reaction and its kinetics was driven by the stochastic simulation algorithm (SSA) [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]. These models do not account for one important aspect of the kinetics of gene expression, namely that it consists, as mentioned, of a sequential process whose intermediate steps, once initiated, take considerable time to be completed (see e.g. [START_REF] Kandhavelu | Singlemolecule dynamics of transcription of the lar promoter[END_REF]). This feature can be accounted for by introducing 'time delays' in the appearance of the products modeling the process [START_REF] Bratsun | Delay-induced stochastic oscillations in gene regulation[END_REF][START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF][START_REF] Ribeiro | A general modeling strategy for gene regulatory networks with stochastic dynamics[END_REF]. To cope with this necessity, biochemical reactions with stochastic delays were introduced, represented in the following form: R :

n i React i k -→ m j P rod nd j + m h P rod d h (dist h )
where React i , P rod nd j and P rod d h denote respectively the i-th reactant, the jth non-delayed product and the h-th delayed product (n i , m j , m h ∈ N being the stoichiometric coefficients) and dist h denotes the distribution for the delayed introduction of h-th delayed product. Thus, for example, a reaction like A + B k -→ A + C(δ(τ )) represents an interaction between a molecule of A and one of B, that results in the transformation of the reacting B molecule into a C molecule (while the A molecule is left unchanged). The beginning of such interaction is assumed to occur after an exponentially distributed delay proportional to kinetic rate constant k. Once A and B interaction begins, the reacting A molecule is immediately released whereas the produced C molecule will appear after τ time units (i.e. δ(t) denotes a deterministic distribution with duration t). To deal with delayed reactions different adaptations of Gillespie's Stochastic Simulation Algorithm (referred to as "delayed SSA") have been introduced (e.g. [START_REF] Bratsun | Delay-induced stochastic oscillations in gene regulation[END_REF][START_REF] Barrio | Oscillatory regulation of hes1: Discrete stochastic delay modelling and simulation[END_REF][START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF][START_REF] Bratsun | Delay-induced stochastic oscillations in gene regulation[END_REF][START_REF] Ribeiro | A general modeling strategy for gene regulatory networks with stochastic dynamics[END_REF]) yielding to the realisation of dedicated software tools such as SGNSim [START_REF] Ribeiro | SGN Sim, a stochastic genetic networks simulator[END_REF].

Single gene expression model

We consider a model of single gene expression that follows the approach proposed in [START_REF] Ribeiro | A general modeling strategy for gene regulatory networks with stochastic dynamics[END_REF]. Our model differs in that transcription is modeled as a 2step process so as to accurately account for the open complex formation and promoter escape [START_REF] Kandhavelu | Singlemolecule dynamics of transcription of the lar promoter[END_REF]. Each of these processes duration follows an exponential distribution, hence each reaction occurs with an exponentially distributed delay proportional to her kinetic rate (i.e. k t for R 1 , λ 1 for R 2 , etc.). The gene expression system we refer to consists of the following reactions9 :

R1 : P ro + * RNAp k t -→ P rox (6)
R2 : P rox

λ 1 -→ P ro + RBS + RNA(Γ(G len , 0.09)) (7)
R3 : * Rib + RBS k tr -→ RBS(δ(τ1)) + Rib(Γ(G len , 0.06)) + P (Γ(τ5 sh , τ5 sc )) ( 8)

R4 : RBS rbsd -→ ∅ (9)
R5 :

P ro + Rep kr -→ P roRep (10)
R6 : P roRep kunr -→ P ro + Rep (11)
Reactions ( 6) and ( 7) model transcription. In ( 6), an RNAp binds to a promoter (Pro), which remains unavailable for more reactions until reaction (7) occurs. Following reaction [START_REF] Baier | Model checking action-and state-labelled Markov chains[END_REF], which models the promoter escape, both the promoter and the RBS become immediately available at completion of R 2 whereas a complete RNA is released in the system after a further delay established according to a Gamma distribution Γ(G len , 0.09) whose parameter G len is determined by the length of the gene, here set to 1000 nucleotides, and the time spent by the RNAp at each nucleotide, which follows an exponential distribution with a mean of 0.09s [START_REF] Makela | Stochastic sequence-level model of coupled transcription and translation in prokaryotes[END_REF]. Furthermore the kinetic rates k t and λ 1 of reactions ( 6) and ( 7) are set to 1/400s, following measurements for the lar promoter in charge of transcription in Escherichia Coli [START_REF] Kandhavelu | Singlemolecule dynamics of transcription of the lar promoter[END_REF]. Notice that, in this model, RNA will not be substrate to any reaction, and is only modeled as a means to count the number of RNA molecules produced over a certain period of time. In our model, according to the delayed SSA, the time necessary for any reaction to occur follows an exponential distribution whose mean is determined by the product between the rate constant of the reaction with the number of each of the reacting molecules present in the system at that moment.

In Prokaryotes, translation can begin as soon as the ribosome binding site (RBS) region of the RNA is completed. In reaction [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF], a ribosome (Rib) binds to the RBS and translates the RNA. The RBS becomes available for more reactions after τ 1 s. The ribosome is released after Γ(G len , 0.06) seconds. The initiation rate, k tr is set to 0.00042 s -1 [START_REF] Zhu | Studying genetic regulatory networks at the molecular level: Delayed reaction stochastic models[END_REF]. Following measurements from E. coli, we have set τ 3 = 2 s, and Γ(G len , 0.06) to follow a gamma distribution dependent on the gene's length, where each codon is added following an exponential distribution with a mean of 0.06 s [START_REF] Makela | Stochastic sequence-level model of coupled transcription and translation in prokaryotes[END_REF]. Finally, Γ(τ 5 sh , τ 5sc ) is such that it accounts for the time that translation elongation takes, as well as the time it takes for a protein to fold and become active. In this case, we used the parameter values measured from GFP mutants commonly used to measure gene expression in E. coli [START_REF] Megerle | Timing and Dynamics of Single Cell Gene Expression in the Arabinose Utilization System[END_REF]. Finally, we consider also three additional reactions representing respectively: RBS decay (9), promoter repression [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] and its reverse [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF]. Initially, the system has 1 promoter and 100 ribosomes. In the remainder of the paper we illustrate a thorough analysis of the above described single gene model by means of HASL model checking. The single-gene model described by equations [START_REF] Aziz | Model-checking CTMCs[END_REF] to [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] is encoded as an NMSPN by the net depicted in Figure 18. The net includes a place for each species of the model (i.e. Pro, RNAp, Prox, RNA, RBS, Rib, P, Rep and ProRep) plus a number of auxiliary places representing intermediate stages of delayed reactions (i.e. RNA, RBS, P, Rib). Initial marking of the net is set by means of parameters i rep, i rnap, i rib, which correspond to the chosen initial population of the model (note that the promoter place Pro is initialized with one token, as each gene has one promoter region).

Reactions R1 to R6 of the single gene model correspond to subnets (enclosed in dashed rectangles) in Figure 18. Each such subnet contains either a single exponentially-distributed transition (in case of reactions with non-delayed products i.e. R1, R4, R5, R6) or a combination of exponential and non-exponential transitions (in case of reactions with delayed products, i.e. R2 and R3). For example, subnet R3 in Figure 18 represents the encoding of the translation reaction. It consists of: the translation-start event (i.e. exponentially distributed transition labeled s transl, the RBS release event (i.e. deterministically distributed transition clearRBS ), the ribosome releasing event (gamma distributed transition relRib) and the protein production event (i.e. gamma distributed transition prodP ). Observe that the effect of repressed gene-expression can be promptly analyzed by setting up the repressor's initial population (parameter i rep= l r ): unrepressed configurations corresponds to l r = 0, whereas l r > 0 settings correspond to repressed model where the level of repression is proportional to l r > 0.

Properties of single gene model

Table 3 depicts an excerpt of (informally stated) relevant measures of the single-gene model. They are grouped according to different aspects of geneexpression performance. The corresponding HASL encoding is given in Table 4. We briefly illustrate the automata of Table 4 and the associated HASL expressions: A 1 is designed for measures concerning the occurrences of transc and transl events. It accepts all paths of duration T and uses variables, n 1 and n 2 to maintain the number of transc and transl transitions occurring along a path. Different measures can be assessed through different HASL expressions referred to A 1 including:

φ 1a = (A 1 , E[last(n 1 )]); φ 1 b = (A 1 , E[last(n 2 )]
) and φ 4 = (A 1 , E[last(n 2 )/last(n 1 )]) (see Table 3) and their PDF and CDF equivalent φ 2a to φ 3 b . A 2 accepts only paths along which the number of transcriptions is n 1 = C (within time T ). The probability to have precisely C transcriptions within time T can be computed through specification φ 2a = (A 2 , P ). A 3 is designed for measures concerning the amount of time a gene is repressed. Apart from the usual global clock t it uses a timer t r registering the time gene is repressed, hence it grows ( ṫr = 1) only in location l 1 (i.e. repression is ON, corresponding to a marking of place P roRep > 0), while it remains constant ( ṫr = 0) in location l 0 (i.e. marking of place P roRep = 0). Also note that both l 0 and l 1 are initial locations, which is perfectly legal as their constraints make them mutually exclusive (this way A 3 can be used to analyze both repressed and unrepressed configurations of the model). A 4 measures "how likely it is that, within a transcription interval (i.e. the interval between two occurrences of the transc event), at least N translations have been completed". It uses variables n 1 and n 2 (as above) and n 3 to count how many transcription intervals (along a path) contain n 2 ≥ N translations. The result is stored in p 1 = n 3 /n 1 on acceptance. Note that, in this case, we consider an event-bounded observation window consisting of n1 = N 1 transcription events. Measure φ 5 of Table 3 in HASL terms is φ 5 = (A 4 , E[last(p 1 )]). A 5 is designed for measures of sustainment of translation activity under repression (i.e. φ 7 in Table 3). It uses the following variables: n o counting the number of repression intervals (interval between two repression events) in which translation arrested; t o : measuring the translation time-to-arrest in a repression interval (given that translation arrested); T o timer measuring the cumulated t o . Note that translation arrest corresponds to the absence of tokens in all translation related places of the NMSPN model (Figure 18), corresponding to condition: (RBS = 0∧ RBS = 0∧ P = 0∧ Rib = 0). Locations l 0 , l 1 and l 2 are then associated to the following state conditions of the model: repression is off (l 0 ), repression is on and translation off (l 1 ) and repression is on and translation ongoing (l 2 ). All paths of duration t = T are accepted and the target measure10 is be obtained through expression Z = E[last(T o )/last(n o )].

Experiments.

We assessed the previously described HASL measures through experiments executed with the COSMOS model-checker. For time-bounded measures we have considered (following [START_REF] Ribeiro | SGN Sim, a stochastic genetic networks simulator[END_REF]) T = 2 •10 5 as time horizon which roughly corresponds to 60 cell cycles, considering an average cell cycle of about 55 minutes (i.e. 3300s) in the case of E. coli. All experiments have been run with the following settings concerning confidence interval estimation: confidence-level: 99.99%; interval-width: 0.01. Observe that the throughput of translation is roughly twice as much as that of transcriptions, (both in unrepressed condition, as well as, in presence of repression). This is due to the rates of RNA degradation and translation initiation.

Discussion

We discuss here the contribution, in terms of the HASL based analysis of the gene-expression mechanism, considered in this section.

Some of the properties we have introduced here (Table 3) are just examples of typical performance measures (e.g. specification φ 1 , corresponding to the average number of completed transcriptions/translations and its variants φ 2 and φ 3 which refer to the PDF, respectively CDF, of the average number of completed transcription/translation). φ 1 admits for an equivalent characterization through other stochastic logics, for example through the reward-enriched version of CSL [START_REF] Kwiatkowska | Stochastic model checking[END_REF] (supported by PRISM) which provide the user with dedicated reward-operators (e.g. R =? ) and allows one to assess the average value of state and/or transition rewards cumulated along paths. On the other hand φ 2 and φ 3 refer to the PDF, respectively CDF, of a reward-based quantity (i.e. the average number of completed transcription/translation). The reward-based extension of PRISM CSL seems not to admit for a direct equivalent expression of measures like φ 2 , φ 3 as in its syntax the argument ϕ of a probabilistic formula P = ?(ϕ) can only be a time-bounded path-formula (e.g. ϕ ≡ (Φ U t∈I Ψ) where I ⊆ R ≥0 is the time-bounding interval) hence is not affected by reward-structures. On the contrary a logic like CSRL [START_REF] Cloth | Model checking Markov reward models with impulse rewards[END_REF], which also extend CSL reasoning to CTMC enriched with state/impulse rewards, may admit an equivalent formulation for assessing the CDF of a reward-quantity (specifications like φ 3 ) as the argument ϕ of a CSRL probabilistic formula P = ?(ϕ) is a time/reward-bounded path formula (i.e .(Φ U t∈I r∈J Ψ)). The remaining specifications φ 4 , φ 5 , φ 6 , φ 7 (Table 3) are instead examples of properties that cannot be directly encoded through logics based on "simple" state/transition reward-structures like, for example, CSRL and the rewardextension of CSL featured by PRISM. This is because all of φ 4 , φ 5 , φ 6 , φ 7 are based on rewards whose value must be updated in function of conditions observed along a path (e.g. with φ 4 we need to count the number of translation events occurred between two consecutive transcription events) and not only associated to the current state or to the transition being observed as it is the case with the reward-structures supported by CSRL and CSL.

The characterization of path-dependent reward-functions is normally achieved through association of an automaton (which selects the paths with relevant dynamics) to the input (CTMC) model. Thus, in principle, properties such as φ 4 , φ 5 , φ 6 , φ 7 could still admit an equivalent formulation through logics such as CSRL and CSL at the cost of having, for each such property, to manually encode and incorporate an automaton into the considered input CTMC/MRM model.

To avoid such a costly burden, one can resort to automata-based logics (e.g. [START_REF] Donatelli | Model checking timed and stochastic properties with CSL T A[END_REF][START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF][START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF][START_REF] David | Statistical model checking for networks of priced timed automata[END_REF]), i.e. logics that use a formalism using automata and allowing the user to explicitly characterize the dynamic behavior to be analyzed through an automaton. From the point of view of the expressiveness these logics are similar. However, differences exist depending on the class of automata a logic is based on. For example with timed-automata (TA) logics, such as CSL TA [START_REF] Donatelli | Model checking timed and stochastic properties with CSL T A[END_REF] or its n-clock extension [START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF], one can assess the probability that the considered model exhibits sophisticated timing conditions. However, in their original formulation these logics do not support reward-based analysis hence they are constrained to the evaluation of probability measures. Thus, even if semantically the paths associated to timing-dependent properties such as φ 6 ≡ (A 3 , E[(100/T ) * last(t r ))]) and φ 7 ≡ (A 5 , E[last(T o )/last(n o )] could be characterized through CSL TA and its [START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF] extension, the actual target measures cannot since both E[(100/T ) * last(t r ))] and E[last(T o )/last(n o )] are in fact reward measures. Moreover, the target measure E[last(T o )/last(n o )] depends of an integer variable (n 0 , the number of observed transcription intervals in which translation has completed), hence could not be expressed anyway by any logic based on a "pure" TA formalism as TA are limited to clock variables (i.e. variables that can increase with gradient 1 or stay constant and whose update can only be a reset).

Similar considerations apply also to specifications φ 4 and φ 5 which represent a reward-measure (φ 4 ), respectively a probability measure (φ 5 ) based on counting certain events within an observed event pattern. On the other hand, logics based on generalized version of TA, such as Linear Hybrid Automata [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] or Priced Timed Automata [START_REF] David | Statistical model checking for networks of priced timed automata[END_REF], avoid the limitations of TA logics by allowing the modeler to use/combine any type of variables (real, integer, boolean) in order to express the relevant characteristic of the paths to observe as well as the target measure (probability of reward) to be evaluated.

Finally, we stress that the increased expressiveness achieved through automatabased specifications comes at the cost of having to get acquainted with expressing a measure through an automata. This is the main drawback of automatabased specification formalisms. Within HASL the development of a higher level specification language (e.g. based on some sort of temporal-logic-like syntax), suitable for automatic translation into a corresponding automata specification is certainly conceivable and part of the future work, as discussed in [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF].

Conclusion

In this paper we have considered the application of an expressive stochastic specification language, namely HASL, to the analysis of biological models. First we have considered the problem of analyzing oscillations, and in that respect we have presented HASL specifications for assessing both noisy periodicity and noisy alternance. We have applied this approach on a real biological oscillator, evaluating expected values of interesting measures, as well as their PDF and CDF. In the second part of the paper we focused on another type of problem, namely the analysis of gene-expression with stochastic delayed dynamics. This type of models often involve non-exponential delay distribution, and thus cannot be analyzed using numerical stochastic model checking approaches. In the analysis of gene-expression, we demonstrated the potential of HASL by illustrating examples of sophisticated analysis, such as the measurement of the efficiency of the transcription-phase with respect to the translation-phase of the expression process (i.e. by means of measures of the sustainment of translation with respect to transition).

16.67 38.89

Table 5: Configurations of the square wave oscillator: regular waves correspond to a single maximal priority (and in this case the weights are irrelevant); irregular waves correspond to at least two maximal priorities (and in this case the period duration is probabilistic and depends on the corresponding weights). The (exact) mean value (µt p ) and fluctuation (ft p ) of the oscillation period computed for d = 5 are indicated in the last two columns.

the deterministic-timed transitions (represented as thick bars), each of which is associated with a specific delay-duration (i.e. a multiple of the parameter d).

Finally, the possibility of switching between a regular square wave configuration and an irregular one is achieved by means of the parameters of the 3 immediate transitions t 1 , t 2 and t 3 (represented as thin bars) which represent the selection of the duration of the high-level interval. Each such transition is associated with a priority and a weight ((π i , w i ), i ∈ {1, 2, 3}) parameter. A transition with maximal priority is always chosen first. As a consequence a regular square wave is obtained by any (deterministic) configuration which guarantees that only one amongst t 1 , t 2 and t 3 can be chosen (i.e. this is achieved by having a single maximum priority π i > max(π j , π k ) i, j, k ∈ {1, 2, 3} and i = j = k).

Thus we have three possible regular square waves: reg 1 : (π 1 , π 2 , π 3 ) = (1, 0, 0), reg 2 : (π 1 , π 2 , π 3 ) = (0, 1, 0) and reg 3 : (π 1 , π 2 , π 3 ) = (0, 0, 1) On the other hand, irregular square waves correspond to configurations where at least two transitions have equal maximum priority (i.e. π i = π j > π k or π i = π j = π k ). In this case the maximally prioritized transitions are probabilistically selected according to their respective weights. For example with priorities (π 1 , π 2 , π 3 ) = (1, 1, 1) transition t i is selected with probability w i / w j . There are four possible irregular waves listed in Table 5. The mean value µ tp and the fluctuation of the period f tp for an irregular configuration can be straightfor-wardly expressed analytically in function of the priority and weight parameters of transitions t 1 , t 2 , t 3 . For example for the irr 4 : (π 1 , π 2 , π 3 ) = (1, 1, 1) configuration they are given by [START_REF] Ballarini | Query-based verification of qualitative trends and oscillations in biochemical systems[END_REF], respectively [START_REF] Ballarini | Expressive statistical model checking of genetic networks with delayed stochastic dynamics[END_REF].

µ
In the following experiments, we consider the configurations illustrated in Table 5. Note that, in the table, irregular configurations are defined so that all possible durations of the period are equally likely (i.e. the weights of maximally prioritized transitions have all the same value). This induces a (specific) discrete distribution of probability over the oscillation period t p , example with configuration irregular1 (i.e. priorities (1, 1, 0) weights: (1, 1, 0)): P r(t p = 10) = P r(t p = 15) = 0.5. As a consequence the exact mean value (µ tp ) and fluctuation (f tp ) of the period t p can be straightforwardly obtained analytically: the actual values are indicated in Table 5. As expected, deterministic configurations have no fluctuation. On the other hand the most irregular waves (i.e. those whose period varies the most) correspond to priorities configuration (1, 0, 1): such waves will alternate between equally likely high-level phases of duration d and others of duration 4d. Experiments on the square wave oscillator. Figure 24 and Figure 25 report on measurements of the average period, respectively the fluctuation of the period, for the square wave oscillators, obtained through the LHA for measuring noisyperiodic traces (i.e. Figure 7). For these experiments, knowing a priori the amplitude of the square waves (i.e. min=1, max=10), we manually set the LHA thresholds to L = 2 and H = 9. The goal was then to assess measures related to the period. For the null-fluctuation (regular) configurations of the square wave oscillator (i.e. reg1, reg2 and reg3), the measured average duration and fluctuation of the period precisely match the theoretical ones depicted in Table 5. Figure 24 and 25, instead, refer to the irregular configurations. In particular the average value (respectively the fluctuation) of the period is measured as a function of either weight w 2 or w 3 (depending on the considered configuration), hence as a function of the probability of selecting either 2d or 4d as the duration of the high-level phase. For example the curves for configuration irr 2 , irr 3 and irr 4 refer to measuring the effect on the average value and fluctuation of the period of the variation of w 3 ∈ [1, 50] while keeping the w 1 and w 2 fixed. The curves for configuration irr 1 refer to measuring the effect of the variation of w 2 ∈ [1, 50] while keeping w 1 = 1. For irr 2 , we observe that augmenting w 3 increases the probability of observing a period of length 5d, leading the oscillator to behave, for w 3 → ∞, like the regular configuration reg 3 . This is confirmed by the fact that the irr 2 curves in Figure 24 and 25 tends towards the µ tp = 25, respectively s 2 tp = 0 values of the regular configuration reg 2 . Notice that for all plots in Figure 24 and 25 the measured value on x = 1 matches the exact value of the corresponding irregular configuration (see Table 5), which is a good indication of the accuracy of the HASL measure. For the measure period fluctuation Figure 24 notice that the decreasing shape of each curve confirm the intuition: the irregularity of the period (i.e. its fluctuation) depends on the distribution of probability between the possible different durations. Intuitively the more uniform such probability is distributed amongst several possible durations the higher the fluctuation of the period. By contrast at the limit w 2 → ∞ (or w 3 → ∞) the fluctuation asymptotically tends to zero. In a similar fashion, we also validated the LHA for assessing noisy-periodically alternating traces. In particular with the LHA in Figure 5 we measured the amplitude and the period of the square wave oscillator obtaining the correct results (to save space we omit to report them here).

A.2. Experiments for measuring the oscillation of the Repressilator

In the following we provide some information concerning the experiments for the analysis of the repressilator model described in Section 3.2.1. All these experiments have been done with COSMOS and using the default configuration for what concerns the confidence-level (i.e. 0.99) and the interval-width (i.e. 0.001). Table 6 refers to experiments for the measure of the PDF of maximal peaks of oscillations obtained through A peaks (see Figure 11).

It shows various relevant information for each experiment performed through A peaks . Specifically for each experiment Table 6 illustrates: the different values of A peaks 's noise parameter used for the experiment (first column); the number of generated path (second column); the actual width of the confidence interval at the end of the estimation (third column); the actual (total) simulation for completing the experiment (fourth column); the number of jobs executed in parallel11 (fifth column); the build time for compiling the simulator (sixth column). COSMOS employs a code-generation scheme through which a customized simulator is built for every single experiment by compilation of the GSPN model and of the LHA specification (see [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF]). The total run time of an experiment is then given by the sum of the simulation time (fourth column) and the time to build (compile) the simulator for the experiment (sixth column). The interested reader can reproduce the experiments discussed in this paper by retrieving the GSPN models and HASL specifications files available at COSMOS web-page http://www.lsv.ens-cachan.fr/Software/cosmos/. Instructions

Figure 1 :

 1 Figure 1: Example of NMSPN: model of reaction R 1 and R 2 of single-gene system

Figure 2 :

 2 Figure 2: Examples of LHA for selecting paths of the NMSPN model of Figure 1

Example 2 . 3 .

 23 let us consider A a and A b the two LHAs of Figure 2(a), respectively Figure 2(b). We provide below examples of HASL specifications obtained by combining A a and A b with appropriate HASL expressions:

Figure 3 :

 3 Figure 3: Ac: detecting N occurrences of transc events within time T in the NMSPN model of Figure 1

Figure 4 :

 4 Figure 4: Example of noisy alternating trace corresponding to a given noise level (δ).

Figure 5 :

 5 Figure 5: A peaks : an LHA for detecting local maxima/minima (for observed species A) of noisy periodic traces with level of noise δ.

Figure 6 :

 6 Figure 6: Example of noisy periodic trace corresponding to a given (L < H) partition of a DESP state space.

  E,(n<N),; E,(n<N ^top = 0),; E,(n<N),{top := 1} E,(n<N),; ],(n = N ),; ] ,( n = N ) ,; ],(n = N ),; A  L L A  H A H E,(n<N),; E,(n<N),; E,(n<N),; E,(n = 1 ^top = 1), {n++, t := 0, top:= 0} E,(2  n<N ^top = 1), {n++, top:= 0} tp := f (tp, n), ft p := g(ft p , tp, n), tp := 0} E,(0  n  1 ^top = 1), {n++, top:= 0, tp := f ( tp, n), tp := 0} low ṫ Skipping the initial transient, until initT second have passed and A<=L A> L l 0 0 ],t initT, ; E , > , { t := 0 , n := 1 } E, >, ; E, >, ; A per end ],t initT, {t := 0, n := 1}

Figure 7 :

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

Figure 8 :

 8 Figure 8: The cyclic topology of the repressilator circuit.

R1 1 :

 1 P ro1 + RNAp k -→ RNAp + P ro1 + P1 R2 1 :P1 + P ro2 kr -→ P1P ro2 R1 2 :P ro2 + RNAp k -→ RNAp + P ro2 + P2 R2 2 :P2 + P ro3 kr -→ P2P ro3 R1 3 :P ro3 + RNAp k -→ RNAp + P ro3 + P3 R2 3 :P3 + P ro1 kr -→ P3P ro1

Figure 9 :

 9 Figure 9: GSPN model of the three-genes repressilator: places P1, P2 and P3 represent the population of the three expressed proteins.

Figure 11 :

 11 Figure 11: PDF of the local maxima of repressilator's oscillation peaks for different protein's decay rate and for different values of the accepted noise (i.e. δ) parameter of A peaks .

  (a) for example). Here we propose that appropriate values of L and H should be chosen as a function of the distribution of probability of the detected local 25

Figure 12 :

 12 Figure 12: PDF of the local minima of repressilator's oscillation peaks as a function of the decay rate and for different values of the accepted noise.

  depicts the measures of the average height of maximal peaks (Figure14(a)) and of the average period duration (Figure 14(b)) computed through A peaks , respectively A per for d 1 ∈ [0.005, 0.02]. The measured height of maximal peaks (assessed for different values of the noise parameter of A peaks ) monotonically decreases thus confirming the intuition that peaks shrink as P 1 degradation speeds up.

  average duration of a period.

Figure 14 :

 14 Figure 14: Average height of maximal peaks (left) and period duration (right) measured with A peaks , respectively, Aper on the repressilator with initial RNAp = 30 and in function of the degradation rate of protein P 1

Figure 15 :

 15 Figure 15: traces of the repressilator corresponding to a 33x faster protein expression (RNAp =1000) and for different values of protein's degradation rate d (i.e the speed of the degradation of the protein P i when attached to the promoter).

  detail for d ∈ [0.01, 0.25]

Figure 16 :

 16 Figure 16: Average period duration versus period fluctuation for the repressilator with 33x faster expression measured in function of degradation rate d .

Figure 17 :

 17 Figure 17: Convergence of the mean value of the cycle duration and of the maximal/minimal peaks for the repressilator and with respect to the number of observed cycles/peaks cycle duration (Figure 17(a)) and maximal/minimal peaks (Figure 17(b)) in function of N , the number of observed cycles, resp. peaks. Each graph contains different plots corresponding to different settings of the relevant parameters of the experiments (i.e. the L, H thresholds for the cycle measures of Figure 17(a) and the noise level for the peaks measures of Figure 17(b)). The obtained results indicate that, both for the cycles and peaks measures, a stable value is reached at

Figure 18 :

 18 Figure 18: NMSPN model of Single Gene system with delayed stochastic dynamics

Figure 19 :Experiment 1 .Experiment 2 .

 1912 Figure 19: PDF and CDF of completed transcriptions within T

  molecules of repressor percentage of repressed-time v translation-off-time % rep-time % transl-off-time num-transc num-transl (b) % rep-time v transc and transl

Figure 20 :

 20 Figure 20: Exp. transcriptions and translations and percentage of repressed time

Figure 20 (

 20 b), plots two measures of timing: the percentage of time gene is repressed (A 3 ) and the percentage of time no translation activity is going on (variant of A 3 ) when system is observed for duration T and in function of the level of repression (num. of repressor molecules on the x-axis). To observe also the trend of transcription and translation activity in function of repression level, Figure20(b) also includes two curves referred to the expected number of transcriptions, respectively translations, within T . Observe

Figure 21 :

 21 Figure 21: Measure related to translation activity

Figure 23 :

 23 Figure 23: traces of a square wave oscillator with: (left) constant period (∆p = 2d); (right) probabilistically variable period (∆p ∈ {2d, 3d, 5d})

Figure 24 :

 24 Figure 24: Measured average duration of the period for different configurations of the square wave oscillator, as a function of the weights, i.e. for irr 1 w 2 ∈ [1, 50] and w 1 = 1, for irr 2 w 3 ∈ [1, 50] and w 2 = 1, for irr 3 w 3 ∈ [1, 50] and w 1 = 1 and for irr 4 w 3 ∈ [1, 50] and w 2 = w 1 = 1.

3 Figure 25 :

 325 Figure 25: Measured fluctuation of the period duration for different configurations of the square wave oscillator, as a function of the weights i.e. for irr 1 w 2 ∈ [1, 50], w 1 = 1, for irr 2 w 3 ∈ [1, 50] and w 2 = 1, for irr 3 w 3 ∈ [1, 50] and w 1 = 1 and for irr 4 w 3 ∈ [1, 50] and w 2 = w 1 = 1.

  x < A-δ) or a

			Data variables
	name	domain update definition	description
	t	R ≥0	reset	time elapsed since begin-ning measure (first non-
				spurious period)
	n M (n m )	N	increment	counter of detected lo-
				cal maxima (n M ), minima
				(n m )
	x	N	current value of	(overloaded) variable stor-
			observed species A	ing most recent detected
				maximum/minimum
	Smax(Smin)	N		sum of detected maxima
				(minima)
	Lmax[](Lmin[])

Table 3 :

 3 Properties of the Single Gene model

  that the presence of a single repressor is sufficient for the gene to remain repressed for 83% of the time and, likewise, for translation to be non-existing for 85% of the observation time (whereas in absence of repressor, translation activity is only non-existing for about 4% of the time). Figure21(a)shows the PDFs of random variable n 2 : number of completed transcription within a transcription interval (i.e. within two consecutive transcription completions) corresponding to query φ 4 : (A 4 , E[last(p 1 )]. This is computed for the unrepressed model and for two configurations of the repressed model (rep(1) and (rep(2)). Outcomes indicate that, in presence of repression, the probability density is more "distributed" than the bell shaped one corresponding to the unrepressed configuration. Furthermore, increasing the level of repression seems to have no great impact on the probability density (plots rep(1) and rep(2) are essentially identical).

	Experiment 4.

  [START_REF] Arkin | Stochastic Kinetic Analysis of a Developmental Pathway Bifurcation in Phage-l Escherichia coli[END_REF] , E[last(n of f )/last(n rep )]) in function of the RBS decay rate (rbsd). We conducted our experiment with rbsd in the interval [0.001, 4], which includes rbsd = 0.01 i.e. the value complying with experimental evidence used in the"standard" model's configuration. Obtained results indicate quite sensibly that translation sustainment is inversely proportional to RBS decay. It should be noted that with rbsd < 0.004 the translation sustainment is actually increasing with rbsd (not very evident in plot of Figure21(b)). This is because, by definition, query φ 5 : (A 5 , E[last(n of f )/last(n rep )]) measures the sustainment of translation on condition that sustainment lasts lesser than repression. With rbsd < 0.004, however, decay is so slow that with high probability sustainment lasts longer than repression, while with low probability it lasts less. Thus it is sensible that the duration of translation sustainment not exceeding repression duration increases for very small rbsd, i.e. rbds ∈ [0, 0.004].

Table 6 :

 6 Experiments for measuring the PDF of maximal peaks through A peaks Data concerning experiments run with A peaks for measuring the PDF of the maximal peaks of oscillations for the repressilator model

	noise	tot. paths	interval width	sim. runtime (s)	jobs	build time (s)
	1	23000	0.000224	109.08	4	1.12
	2	15800	0.000326	166.91	4	1.44
	3	15650	0.000330	282.13	4	1.34
	5	15250	0.000339	466.35	4	1.25
	10	15050	0.000344	750.84	4	1.10
	15	16400	0.000316	969.44	4	1.37
	20	18800	0.000302	1388.56	4	1.40
	25	25150	0.000209	2753.07	4	1.15

This is useful for eliminating the effect of the initial transient from long run measures as already discussed in[START_REF] Amparore | Expressing and computing passage time measures of gspn models with hasl[END_REF].

assuming the oscillator starts in low: slight variants of (1) hold in case the oscillator starts in either mid or high.

the choice of beginning measuring in low is arbitrary, equivalently Aper can be defined so that analysis starts in either mid or high.

again equivalent versions of Aper can be easily obtained which detect periods by considering different starting points, e.g periods that start with a crossing from mid to high.

in fact it requires verifying that the steady-state measure for certains states of the timedautomaton expanded CTMC add up to 1.

peaks with height below the chosen noise are not detected by A peaks hence their probability is correctly zero.

note that symbol * prefixing a species name in the above reactions means that the reactant is not consumed in the reaction. This is applied for simplicity to those reactants such as ribosomes, which exists in large amounts, and thus fluctuations in their numbers will not be significant in the propensity of reactions.

note that with a time-bounded measurement, as with A 5 , measuring may stop in any instant (not necessarily at the end) of a repression interval: this is not a problem as To and no are updated only when translation arrests, thus if bound T is reached before translation arrests, measure To/no will correctly refer to the duration of translation sustainment over all completed repression cycles

i.e. with COSMOS the user can choose to split the simulation of batches of traces in parallel on a given number of jobs which will be then mapped by the OS on the actual cores of the CPU.

Acknowledgments

The authors would like to warmly thank Andre S. Ribeiro, Jarno Mäkelä and Ilya Popatov for their precious contribution to this work, particularly with respect to the part regarding modeling of gene-expression with stochastic delayed dynamics.

,(n 1 =N 1 ),

sustainment of translation under repression

Table 4: LHA for various measures of the Single Gene model A.1. Checking the validity of the measures of period obtained through A per on a simple example of probabilistic square wave oscillator

In order to validate the HASL based approach for oscillation related measures, we consider a simple example of probabilistic square wave oscillator, corresponding to the NMSPN model given in Figure 22. The interest in such a toy-example of synthetic oscillator is twofold. First its simplicity makes it possible to estimate measures (period, fluctuations) analytically and thus compare them with the measures obtained with HASL. Then its oscillatory characteristics (i.e. amplitude and period of oscillations) are well-defined and can be easily tuned through the model's parameters thus swapping between a deterministic oscillator (i.e. constant period) and a probabilistic one (period duration's associated to a discrete probability distribution). The model consists of a single biochemical species A whose population (corresponding to the marking of place A in Figure 22) periodically alternates between 1 molecule (low-level) and 10 molecules (high-level) through step-jump transitions. The duration of the low-level phase is fixed (d) and determined by the (deterministically distributed) timed-transition T u , where d is a parameter of the NMSPN. On the other hand, each high-level phase may have any of 3 possible durations: d, 2d or 4d. The different durations available for the highlevel phase are determined through the model's configuration, i.e through the chosen values for the priorities (π 1 , π 2 , π 3 ) and the weights (w 1 , w 2 , w 3 ) of the immediate transitions t 1 , t 2 and t 3 (Table 5). According to the chosen configuration, we may distinguish between regular square waves (e.g. left-hand side in Figure 23), corresponding to waves of constant period (i.e. either ∆p = 2d, ∆p = 3d or ∆p = 5d), as opposed to irregular square waves (e.g. right-hand side in Figure 23), i.e. waves whose period vary probabilistically.

The NMSPN model of the square wave oscillator in Figure 22 consists of 3 major places: A (representing the population of species A) plus the two mutually-exclusive places HIGH and LOW (representing the current level of the wave). The timing of the wave (i.e. the period duration) is driven by on how to install and run COSMOS are available in the archive containing the entire COSMOS distribution (available on the COSMOS web-page).