<figure><img src="image_19.png" /><figcaption>Figure 3.7: Top: bi-dimensional weighted histograms of posterior distribution $\pi_{\varphi_i-ABC}$ for ER system obtained with Automaton-ABC-SMC using constant $k_3 = 1$, prior distributions $\pi_{k_1}(.)$, $\pi_{k_2}(.) \sim U(0, 100)$ and using 1000 particles. Mid: kernel density estimation of the satisfaction probability function. Bottom: validation through estimation of the satisfaction probability function by statistical model checking (95% confidence semi-interval with half-width 0.01).</figcaption></figure>
# Chapter 4

## Formal analysis of stochastic
oscillators

### 4.1 Introduction

Many real life systems are characterised by some form of oscillatory dynamics examples of which can be found in many domains, ranging, from ecology, e.g. with the early Lotka-Volterra *prey-predator* model [Lot09, Vol28], to social systems, where periodic behaviour can be observed in many social context for example with the arrivals of patients at hospital's emergency department [BDHA20] and, specially, in biology where oscillations are prominent dynamics at the core of many fundamental biological processes [Gol02].

Here we focus on two complementary aspects related to modelling oscillators: the assessment of meaningful indicators (e.g. the period and amplitude of oscillations) from an oscillator model and the *reverse engineering* of an oscillator (i.e. how to tune an oscillator model so that it meets a given a target oscillatory behaviour). As for the first aspect, the analysis of mathematical models of oscillators is a well established problem in the literature, which typically boils down to *limit-cycle analysis* for *deterministic oscillators* (i.e., oscillators modelled as a systems of Ordinary Differential Equations), or, alternatively is achieved by application of signal processing methods such as Fast Fourier Transformation (FFT) or autocorrelation analysis. In this contribution we focus on discrete-state stochastic oscillators and propose an alternative approach, which differently from classical signal processing, is based on finding a *temporal logic* characterisation of periodicity for noisy signals issued by a stochastic model. With such a formal characterisation of what we name *noisy periodicity* we can then take advantage of the model checking framework so to automatically assess the period and amplitude of oscillations.

As for the second aspect we look at analysing the relationship between the model's parameters and the character of the oscillations. More specifically we tackle the problem of how to *reverse engineer* an oscillator, that is, given a target oscillatory behavior, e.g.
a desired period of oscillation, we aim at inferring the model’s parameters that *match* the desired behaviour. We obtain such a framework based on the automaton-ABC approach outlined in Chapter 3.

**Outline of the contribution.** This chapter report on the contribution in the matter of formal analysis of stochastic oscillators. We demonstrate the effectiveness of the HASL formalism as a means for specifying and automatically estimate oscillation related measures, such as the period and amplitude of noisy oscillators. We do so by introducing two alternative characterisations of noisy periodicity, one relying on a pair of selected observational thresholds, the other based on the detection of (noisy) local maxima and minima. We then provide a formal encoding in terms of linear hybrid automata (LHA) of the HASL logic for both characterisations of periodicity, which we demonstrate being capable of detecting the periods, respectively the peaks, of the stochastic oscillator they target and to estimate *on-the-fly* “classical” indicators such as the average duration and the average amplitude of the oscillations as well as more sophisticated ones like the period variance, which allow for assessing the *regularity* of the oscillator. Differently from Spieler’s approach [Spi13], which allows one for qualitatively establishing whether a CTMC model is a sustained oscillator and, if so, to assessing some quantitative characteristics (e.g. the average period duration), the methodology we introduce here is limited to assessing quantitative characteristics of a stochastic model which is known (or believed to) oscillates.

## 4.2 Formal analysis of stochastic oscillators

Intuitively an oscillation is the periodic variation of a quantity around a given value. In mathematical terms this is associated with the definition of (non-constant) periodic function. i.e. a function $f : \mathbb{R}^+ \to \mathbb{R}$ for which $\exists t_p \in \mathbb{R}^+$ such that $\forall t \in \mathbb{R}^+$, $f(t) = f(t + t_p)$, where $t_p$ is called the period as in e.g., Figure 4.1(a). In the context of stochastic models such a “deterministic” characterisation of periodicity is of little relevance, as the trajectories of a stochastic oscillator being strictly periodic (as in $f(t) = f(t + t_p)$), will have (unless in degenerative cases) zero probability. More generally the paths of (discrete-state) stochastic oscillators are characterised by a remarkable level of noise: period occurrences differ in duration as well as in amplitude as in e.g., Figure 4.1(b).

The difficulty we face is that of finding a suitable characterisation of what should be considered as a *period realisation* which agrees with the intuitive notion of periodicity, while copying with the noisy nature of stochastic oscillator paths. The relevant question in this respect is “what shall one consider as a period realisation”? We propose two complementary approaches to determine what is a noisy period realisation. With the first, inspired by Spieler [Spi13] and further developed in [BD15], one has to consider a partition of the domain of the oscillating quantity in three region: two “extreme” regions (*low*, respectively, *high*) separated by an intermediate region (*mid*). A period realisation is then considered as the delay between the beginning of two consecutive sojourns in either of the
<figure><img src="image_20.png" /><figcaption>(a) deterministic oscillation: period occurrences and peaks are deterministic (period equal to 2, maxima at 2 and minima at 0)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) noisy oscillation: the period occurrences depend on the “observation point” (blue line) and are noisy.</figcaption></figure>

<caption>Figure 4.1: Deterministic versus stochastic oscillations</caption>

extreme regions, interleaved by a sojourn in the opposite extreme region (see Figure 4.2(b)). With the second one, introduced in [BD15, Bal15a], a period realisation is considered as the delay between two consecutive local maxima (or equivalently local minima), bearing in mind that the characterisation of local maxima and minima in a noisy oscillator is not trivial as one has to be able to distinguish between actual local maxima/minima points and critical points that correspond to noisy spikes.

### 4.2.1 Temporal logic characterisation of oscillations

The possibility to automatically assess the periodic character of a system has attracted the attention of researchers in model checking community in recent times, the basic idea being to identify suitable temporal logic formulae that can be fed in to a model checker to qualitatively establish the existence of periodicity and (possibly) to quantitatively assess relevant indicators (e.g. period duration, amplitude of oscillations).

In an early work Fages, et al. [CRCD+04] tackled the problem w.r.t. transition system models of cell cycle and proposed the following (branching-time) CTL formula to capture the existence of periodicity:

$$EG((\mathbf{p} \Rightarrow EF(\neg \mathbf{p}))) \wedge (\neg \mathbf{p} \Rightarrow EF(\mathbf{p}))) \qquad (4.1)$$

Formula (4.1) essentially states that there should exist an execution of the model such that *infinitely often* if a condition p is satisfied then it will be not satisfied and vice-versa if it is not satisfied it will be satisfied

In a seminal work with Mardare and Mura [BMM09] we considered a similar CTL-like reasoning to query whether a CTMC population model of a simple 3-way synthetic
<figure><img src="image_22.png" /><figcaption>Figure 4.2: Noisy periodic trajectories of stochastic oscillators</figcaption></figure>

biochemical oscillator exhibits sustained oscillations consisting of the following formula:

$$AG(((X_i = k) \Rightarrow EF(X_i \neq k)) \wedge (((X_i \neq k) \Rightarrow EF(X_i = k)))) \quad (4.2)$$

Formula (4.2) represents a stricter version of (4.1) and establishes, by replacing the outermost existential path quantifier E with the stricter universal quantifier A and by instantiating p with ($X_i = k$), that if the system's evolution reaches a state where the i-th species is $X_i = k$ ($k \in \mathbb{N}$ being a constant) then it has to be possible to reach a future state such that ($X_i \neq k$) and vice versa. Through basics equivalences we then derived a probabilistic CSL counterpart of (4.2), that is:

$$P_{\leq 0}[F(((X_i = k) \wedge P_{\leq 0}[F(X_i \neq k)])) \vee (((X_i \neq k) \wedge P_{\leq 0}[F(X_i = k)])))] \quad (4.3)$$

which can be used to rule out states which allow for converging evolutions, that is states for which there's a positive probability to reach a state $X_i = k$ and never leave it<fnref n="1" />.

The problem with formulae for capturing sustained oscillations of population models, such as (4.2) and its CSL counterpart (4.3), is that they account also for those trajectories whose infinite oscillations around k corresponds to stochastic noise spikes (e.g. trajectories such as $\sigma_i : (x_i=k, x_i=k+1)^\omega$) rather than to actual oscillations therefore we proposed a slight variant formula:

$$AG(((X_i = k) \Rightarrow EF(X_i \leq k-n \wedge X_i \geq k+n)) \wedge (((X_i \leq k-n \wedge X_i \geq k+n) \Rightarrow EF(X_i = k)))) \quad (4.4)$$

which is meant to rule out oscillations corresponding to an arbitrary level of noise n.

<footnote n="1">¹Notice that (4.3) relies on the assumption that all trajectories hit the level of exactly k molecules, that is, the maximum change of the molecule number is one.</footnote>
In a follow up work [BG10] we proposed an evolved qualitative characterisation of periodicity by considering notions of paths' monotonicity, specifically by identifying the monotonic increase, respectively decrease, of the i-th species along a path by enriching the CTMC model with dedicated boolean variables named inc_i, respectively dec_i. That resulted in the following probabilistic LTL formula

$$P_{=?}[\text{inc}_i U(X_i = k \wedge (\text{dec}_i U(X_i = j)))] \quad (4.5)$$

which can be used to measure the probability of paths whose initial state $X_i = j$ and that exhibit a one-period prefix of amplitude $k-j$ consisting of monotonic increase up to $X_i = k$ followed by a monotonic decrease back to the initial level $X_i = j$.

$$P_{=?}[\text{inc}_i U(X_i = k \wedge ((k-n \le X_i \le k+n) \wedge (\text{dec}_i U(X_i = j))))] \quad (4.6)$$

To relax the quite strict characterisation given by (4.5) we introduced the variant formula (4.6) that allows the fluctuation to stay at the peak, that is, around a molecule level of $k$ (module noise band $n$), for an unlimited amount of time, before beginning the monotonic descent to $X_i = j$.

Unfortunately all previously discussed attempts to come up with a temporal logical characterisation of oscillations only partially comply with the start goal of this investigation. For example qualitative CSL formulae like (4.3) although will be unsatisfied by damped oscillators, will also be satisfied by any ergodic CTMC model no matter if the model actually oscillates. On the other hand quantitative CSL formulae like (4.5) and (4.6) help us quantifying how likely oscillations of given amplitude are, but they do not actually help to disambiguate between sustained and damped oscillations as they inherently account for a finite number periods.

A major progress in this respect is given by the PhD work of David Spieler which lead to the publication of paper [Spi13]. The main idea of Spieler was to characterise noisy oscillations by partitioning the state space of the oscillating quantity in 3 regions, a high region where the maxima points of the oscillations are supposedly located, a low region where the minima are located and an intermediate mid region which is traversed in between each consecutive pair of minimum and maximum. Given such characterisation Spieler introduced a corresponding period detector deterministic timed automaton (DTA), i.e. a DTA capable of accepting noisy periodic paths corresponding to the considered low-mid-high state-space partition. Spieler then demonstrated that deciding whether a Markov population model (MPM) admits sustainable oscillations corresponds to check whether a specific steady-state condition holds over the corresponding, so-called, period detector Markov population model (PDMPM), i.e. the (DTA × MPM) product process. Undoubtably a major step forward Spieler's approach suffers nonetheless of some drawback: firstly being based on numerical model checking approaches it is affected by state-space explosion, even more so given the size of the product process PDMPM is at least thrice<fnref n="2" />

<footnote n="2">²The period detector DTA consists of 3 states.</footnote>
the size of the MPM alleged oscillator. Secondly it does not naturally comply with the separation of concern modelling philosophy, as it requires the integration of the detector DTA within the model, something which probabilistic model checking tools (e.g. PRISM) that do not support automata based property languages suffer of. Finally the choice of relying on numerical model checking inherently limits the expressiveness of the oscillation related indicators that one may want to consider: Spieler's outline a transient-analysis based algorithm to approximate the PDF/CDF of the oscillation's period duration however extending such an approach so that other indicators such as, e.g., the variance of the period and/or the PDF/CDF of the oscillation's amplitude looks non trivial if feasible at all.

In the remainder we outline our contribution in this context, specifically, we describe how, by taking advantage of the more expressive hybrid automata based HASL model checking approach, we extended Spieler's DTA period detection analysis so that more sophisticated oscillation related indicators can be assessed.

### 4.2.2 Partition based noisy periodicity

In the remainder we consider $n$-dimensional DESP population models, that is, models whose state space is contained in $\mathbb{N}^n$ and we introduce the notion of *noisy periodicity* w.r.t. to one dimension of the model which, for the sake of simplicity, we often refer to as a species $A$ whose periodicity we want to assess the idea behind such characterisation is to establish a partition (w.r.t. to the $A$ dimension) of the model state-space induced by two threshold levels $L, H \in \mathbb{N}$ with $L < H$ and we say that, a path of the model oscillates (w.r.t. $A$) or, equivalently is *noisy periodic*, if the projection of path over the $A$-dimension infinitely often alternate between a traversal of the $[0, L]$ and of the $[H, \infty)$ regions. This is formalised by Definition 4.1.

**Definition 4.1 (noisy periodic trajectory).** For $\mathcal{M}_\theta$ an $n$-dimensional DESP population model, let $dom(\mathcal{M}_\theta)$ be its state space, and $dom_i(\mathcal{M}_\theta)$ the projection of $dom(\mathcal{M}_\theta)$ along the $i^{th}$ dimension $1 \le i \le n$ of $\mathcal{M}_\theta$. Let $L, H \in \mathbb{N}$, $L < H$, be two levels establishing the partition $dom_i(\mathcal{M}_\theta) = low \cup mid \cup high$ with $low = [0, L)$, $mid = [L, H)$ and $high = [H, \infty)$. A path $\sigma \in Path(\mathcal{M}_\theta)$ is said *noisy periodic* w.r.t the $i^{th}$ dimension, and the considered $L, H$ induced partition of $dom_i(\mathcal{M}_\theta)$ if the projection $\sigma_i$ visits the intervals $low, mid$ and $high$ infinitely often.

**H/L-crossing points.** Given a noisy periodic trace $\sigma_A$ we denote $\tau_{j\downarrow} (\tau_{j\uparrow})$, the instant of time when $\sigma_A$ enters for the $j$-th time the low (high) region. $T_{\downarrow} = \cup_j \tau_{j\downarrow}$ (resp. $T_{\uparrow} = \cup_j \tau_{j\uparrow}$) is the set of all low-crossing points (reps. high-crossing points). Observe that $T_{\downarrow}$ and $T_{\uparrow}$ reciprocally induce a partition on each other. Specifically $T_{\downarrow} = \cup_k T_{k\downarrow}$ where $T_{k\downarrow}$ is the subset of $T_{\downarrow}$ containing the $k$-th sequence of contiguous low-crossing points not interleaved by any high-crossing point. Formally $T_{k\downarrow} = \{\tau_{i\downarrow}, \dots, \tau_{(i+h)\downarrow} | \exists k', \tau_{(i-1)\downarrow} < \tau_{k'\uparrow} < \tau_{i\downarrow}, \tau_{(i+h)\downarrow} < \tau_{(k'+1)\uparrow}\}$. Similarly $T_{\uparrow}$ is partitioned $T_{\uparrow} = \cup_k T_{k\uparrow}$ where $T_{k\uparrow}$ is the subset of $T_{\uparrow}$ containing the $k$-th
sequence of contiguous high-crossing points not interleaved by any low-crossing point. For path $\sigma_A$ in Figure 4.2(b) we have that $T_{\downarrow} = T_{1\downarrow} \cup T_{2\downarrow} \cup T_{3\downarrow} \dots$ with $T_{1\downarrow} = \{\tau_{1\downarrow}, \tau_{2\downarrow}\}$, $T_{2\downarrow} = \{\tau_{3\downarrow}\}$, $T_{3\downarrow} = \{\tau_{4\downarrow}\}$, while $T_{\uparrow} = T_{1\uparrow} \cup T_{2\uparrow} \cup T_{3\uparrow} \dots$ with $T_{1\uparrow} = \{\tau_{1\uparrow}\}$, $T_{2\uparrow} = \{\tau_{2\uparrow}\}$, $T_{3\uparrow} = \{\tau_{3\uparrow}\}$. Based on H/L crossing points we formalise the notion of *period realisation* for a noisy period path.

**Definition 4.2** ($k^{th}$noisy period realisation). For $\sigma_A$ a noisy periodic trajectory with crossing point times $T_{\downarrow} = \cup_{k \ge 1} T_{k\downarrow}$, respectively $T_{\uparrow} = \cup_{k \ge 1} T_{k\uparrow}$, the realisation of the $k^{th}$ noisy period, denoted $t_{p_k}$, is defined as $t_{p_k} = \min(T_{(k+1)\downarrow}) - \min(T_{k\downarrow})^3$.

Figure 4.2(b) shows an example of *period realisations*: the first two period realisations, denoted $p1$ and $p2$, are delimited by the *mid-to-low* crossing points corresponding to the first entering of the *low* region which follows a previous sojourn in the *high* region and their duration (as per Definition 4.2) is $t_{p1} = \tau_{3\downarrow} - \tau_{1\downarrow}$ respectively $t_{p2} = \tau_{4\downarrow} - \tau_{3\downarrow}$. Notice that the time interval denoted as $p0$ does not represent a complete period realisation as there's no guarantee that $t=0$ corresponds with the actual entering into the *low* region. Definition 4.2 correctly does not account for the first *spurious* period $p0$. Relying on the notion of period realisation we characterise the *period average* and *period variance* of a noisy periodic trace. Observe that the *period variance* allows us to analyse the regularity of the observed oscillator, that is, a “regular” (“irregular”) oscillator is one whose traces exhibit little (large) period variance.

**Definition 4.3** (*period mean*). For $\sigma_A$ a noisy periodic trajectory the period average of the first $n \in \mathbb{N}$ period realisations, denoted $\bar{t}_p(n)$, is defined as $\bar{t}_p(n) = \frac{1}{n} \sum_{k=1}^n t_{p_k}$, where $t_{p_k}$ is the $k$-th period realisation.

Observe that for a sustained oscillator, the average value of the noisy-period, in the long run, corresponds to the limit $\bar{t}_p = \lim_{n \to \infty} \bar{t}_p(n)$.

**Definition 4.4** (*period variance*). For $\sigma_A$ a noisy periodic trajectory the period variance of the first $n \in \mathbb{N}$ period realisations, denoted $s_{t_p}^2(n)$, is defined as $s_{t_p}^2(n) = \frac{1}{n} \sum_{k=1}^n (t_{p_k} - \bar{t}_p(n))^2$, where $t_{p_k}$ is the $k$-th period realisation and $\bar{t}_p(n)$ is the period average for the first $n$ period realisations.

We introduce an HASL specification which consists of an LHA automaton, named $A_{per}$, which allows for estimating *on-the-fly*<fnref n="4" /> the period mean and variance, that is, as the noisy periodic traces are generated and scanned by $A_{per}$.

<footnote n="3">³$t_{p_k}$ could alternatively be defined as $t_{p_k} = \min(T_{(k+1)\uparrow}) - \min(T_{k\uparrow})$, that is, w.r.t. crossing into the *high* region, rather than into the *low* region. It is straightforward to show that both definitions are semantically equivalent, i.e., the average value of $t_{p_k}$ measured along a trace is equivalent with both definitions.</footnote>

<footnote n="4">⁴Based on an adaptation of the so-called *online algorithm* [Knu97] for computing the mean and variance out of a sample of observations.</footnote>
<figure><img src="image_23.png" /><figcaption>Figure 4.3: The $\mathcal{A}_{per}$ LHA for selecting noisy periodic traces (with respect to an observed species A) related to partition low = (-∞, L], mid = (L, H) and high = [H, +∞).</figcaption></figure>

**The A<sub>per</sub> automaton.** We introduce a LHA named $\mathcal{A}_{per}$ (Figure 4.3) for assessing the mean value (as per Definition 4.3) and the variance (as per Definition 4.4) of the first $N$ noisy period realisations (as per Definition 4.2) measured on paths of a DESP oscillator $M_{\theta}$. The automaton consists of three main locations **low**, **mid** and **high** corresponding to the regions of the partition of $A$'s domain induced by thresholds $L < H$: location **low** corresponds to region $low = (-\infty, L]$, location **mid** to region $mid = (L, H)$ and location **high** to region $high = [H, +\infty)$. The functioning of $\mathcal{A}_{per}$ is as follows: processing starts in either of the 3 initial locations (**low**,**mid**,**high**) depending on the initial state of the oscillator (variable $n_A$ registers the population of oscillating species $A$ and is initialised through autonomous transitions before unfolding of a path begins). While a path of the considered oscillator is simulated the $\mathcal{A}_{per}$ automaton follows the oscillation by moving between locations **low** and **high**, passing through **mid**, following the profile of the observed species $A$. The completion of a loop from **low** to **high** and back to **low** corresponds to detection of a period realisation (as of Definition 4.2). During period detection a number of relevant information is stored in the data variables of $\mathcal{A}_{per}$ (Table 4.1) which are then exploited for estimating the considered target expressions. The analysis of the simulated trajectory ends by entering location **end** as soon as the $N$-th period has been detected.

**An execution of $\mathcal{A}_{per}$**. Let us assume that initially $\mathcal{A}_{per}$ is in location **low**. From **low** the automaton follows the profile exhibited by the observed population $A$, thus moving to **mid** (and possibly back) as soon as the population of $A$ grows and a state of the $mid = (L, H)$ region is entered (i.e., corresponding to the $L < A < H$ invariant of **mid** location becoming satisfied), and then to **high** (and possibly back) as soon as the population of $A$ enters the $high = [H, +\infty)$ region (corresponding to the $A \ge H$ invariant of **high** location). On
<table><caption>Table 4.1: Variables of the $A_{per}$ automaton.</caption><thead><tr><th>name</th><th>domain</th><th>update definition</th><th>description</th></tr></thead><tbody><tr><td>t</td><td>R≥0</td><td>reset</td><td>total time elapsed</td></tr><tr><td>n</td><td>N</td><td>increment</td><td>counter of detected periods</td></tr><tr><td>t<sub>p</sub></td><td>R≥0</td><td>reset</td><td>duration last period</td></tr><tr><td>t̄<sub>p</sub></td><td>R≥0</td><td>f(t̄<sub>p</sub>, t<sub>p</sub>, n) = <sup>t<sub>p</sub>n·n+t<sub>p</sub></sup>&frasl;<sub>n+1</sub></td><td>period mean</td></tr><tr><td>s<sup>2</sup><sub>t<sub>p</sub></sub></td><td>R≥0</td><td>g(s<sup>2</sup><sub>t<sub>p</sub></sub>, t̄<sub>p</sub>, t<sub>p</sub>, n) = <sup>n-1</sup>&frasl;<sub>n-2</sub>·s<sup>2</sup><sub>t<sub>p</sub></sub> + (<sup>t̄<sub>p</sub>-t<sub>p</sub></sup>&frasl;<sub>n</sub>)<sup>2</sup></td><td>period variance</td></tr></tbody></table>

entering the **high** location the boolean variable *top* is set to true (i.e., *top*=1). This allows then for distinguishing between the **mid-to-low** transitions of kind **mid** $\xrightarrow{E, (...\wedge top=1),...}$ **low**, which correspond to an actual closure of a period realisation (i.e., those $\tau_{j\downarrow}$ preceded by a sojourn in the *high* = [*H*, +∞) region), from those of kind **mid** $\xrightarrow{E, (...\wedge top=0),...}$ **low** which correspond to a return to **low** without having previously sojourned in **high**. Observe that from **mid** location there are four possible (mutually exclusive) ways of entering the **low** location. If the sojourn in **mid** has not been preceded by a sojourn in **high** edge **mid** $\xrightarrow{E, (n<N \wedge top=0),...}$ **low** is enabled. On the other hand if the sojourn in **mid** has been preceded by a sojourn in **high** but **low** is going to be re-entered for the first time (i.e., $n = -1$) then the timer *t* is reset (representing the start time of actual period detection) and the counter of detected periods *n* is set to zero (again representing the actual beginning of counting of period detection). Furthermore if the sojourn in **mid** has been preceded by a sojourn in **high** and the period to be detected is the first one (i.e., $0 \le n \le 1 \wedge top=1$) then we increment the counter *n* of detected period, we reset the flag *top* and update the value of the average duration of detected period $\hat{t}_p$ while we do not update the variable $s_{t_p}^2$ as in order to update the value of the variance of the detected period duration we need that at least two periods have been detected. Finally if the period to be detected is the *n*-th with $n \ge 2$ (i.e., corresponding to guard $2 \le n \le N \wedge top=1$) we do the same update operations of the previous case but also update $s_{t_p}^2$.

The automata uses variable *n* to count the number of noisy periods detected along a trajectory, and stops as soon as the $N^{th}$ period is detected (i.e. event bounded measure). The boolean variable *top*, which is set to true on entering of the *high* location, allows for detecting the completion of a period (i.e. crossing from *mid* to *low* when *top* is true). Two clock variables, *t* and $t_p$, maintains respectively the total simulation time as of the beginning of the first detected period (*t*) and the duration of the last detected period ($t_p$). Finally variable $\bar{t}_p$ maintain the average duration of all (so far) detected periods while $s_{t_p}^2$ stores the variance (or variability) of duration (i.e. how far the duration of each detected period is distant from its average value computed along a trajectory) of all (so far) detected periods.

**Theorem 2.** If a trace $\sigma_A$ is noisy periodic w.r.t. amplitude levels $L, H \in \mathbb{N}$ then it is
accepted by automaton $\mathcal{A}_{per}$ with parameters $L, H$ and $N \in \mathbb{N}$.

*Proof.* See [Bal15b].

<figure><img src="image_24.png" /></figure>

**HASL expressions associated to** $\mathcal{A}_{per}$. We define different HASL expressions to be associated to automaton $\mathcal{A}_{per}$.

• $Z_1 \equiv E[\text{last}(\bar{t}_p)]$: corresponding to the mean value of the period duration for the first $N$ detected periods.

• $Z_2 \equiv \text{PDF}(\bar{t}_p, s, l, h)$: corresponding to the PDF of the average period duration over the first $N$ detected periods, where $[l, h]$ represents the considered support of the estimated PDF, and $[l, h]$ is discretized into uniform subintervals of width $s$

• $Z_3 \equiv E[\text{last}(s_{t_p}^2)]$: corresponding to the variance of the period duration.

Expression $Z_1$ represents the expected value assumed by variable $\bar{t}_p$, that is, the average duration of the first $N$ periods detected along a trace, at the end of accepted trajectory (i.e., a trajectory that contains $N$ periods). Similarly expression $Z_2$ evaluates the PDF of the average duration of the first $N$ periods by assuming the interval $[l, h]$ as the support of the PDF and considering that $[l, h]$ is discretised in $(h-l)/s$ uniform subintervals of width $s$. On the other hand $Z_3$ is concerned with assessing the expected value that variable $s_{t_p}^2$ has at the end of a trace consisting of $N$ noisy periods. By definition (see Table 4.1) $s_{t_p}^2$ corresponds to the *variance* of the duration of the detected periods, (i.e., how much the $N$ periods detected along a trace differ from their average duration). Observe that the measured period variance (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.

### 4.2.3 Peaks-detection based noisy periodicity

The drawback of the partition based period detection scheme discussed in the previous section is that, the detected periods depend on the chosen $L, H$ thresholds, and these have to be chosen by the modeller manually beforehand, i.e., normally by looking at the shape of a sampled trajectory and then choosing where to “reasonably” set the $L$ and $H$ values before executing the measurements with automaton $\mathcal{A}_{per}$. To improve things here we propose a different approach which is aimed at identifying where the peaks (i.e., the local maxima/minima) of oscillatory traces are located.

Since traces of a DESP consist of discrete increments/decrements of at least one unit, it is up to the observer to establish what should be accounted for as a local maximum (minimum) during such detection process. Intuitively a local max/min of a trace $\sigma_A$ is a state $\sigma_A[i]$ ($i \in \mathbb{N}$) that corresponds to a change of trend in the population of $A$. This is formally captured by the following definition.
**Definition 4.5** (local maximum/minimum of a trace). For $\sigma_A$ the projection w.r.t. to species $A$ of a path $\sigma$ of an $n$-dimensional DESP $\mathcal{D}$ population model, state $\sigma_A[i]$ is a maximum, if $\sigma_A[i-1] < \sigma_A[i] > \sigma_A[i+1]$, or a minimum, if $\sigma_A[i-1] > \sigma_A[i] < \sigma_A[i+1]$.

<figure><img src="image_25.png" /><figcaption>(a) local maxima ($M[i]$) and minima ($m[i]$) of
an alternating trajectory $\sigma_A$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) $\delta$-separated local maxima ($M_\delta[i]$) and minima ($m_\delta[i]$) of an alternating trajectory $\sigma_A$.</figcaption></figure>

**Figure 4.4:** Local maxima/minima and $\delta$-separated local maxima/minima.

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

**Definition 4.6** (alternating trajectory). A trajectory $\sigma$ of an $n$-dimensional DESP $\mathcal{D}$ population model is said *alternating* with respect to the $i^{th}$ ($1 \le i \le n$) observed species of $\mathcal{D}$, if $\sigma_i$ contains infinitely many local minima (or equivalently local maxima).

For $\sigma_A$ an alternating trace we denote $\sigma_A^M = M[1], M[2], \dots$, respectively $\sigma_A^m = m[1], m[2], \dots$, the projection of $\sigma_A$ consisting of the local maxima, respectively minima, of $\sigma_A$. Figure 4.4(a) shows the local maxima and minima for an example of alternating trace $\sigma_A$. In the following we point out two simple properties relating the definition of noisy periodic and alternating trace.

**Proposition 5.** If $\sigma_A$ is a *noisy periodic* trace (as of Definition 4.1) then it is also *alternating*.

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

**Proposition 6.** If $\sigma_A$ is an *alternating* trace (as of Definition 4.6) then it is not necessarily *noisy periodic*.

Proposition 6 simply points out that, by definition, an alternating trajectory may be diverging (for example if it consist of increasing steps which are always larger than the decreasing ones), in which case clearly it is not noisy periodic.
In the remainder we introduce a HASL based procedure for detecting the local maxima and local minima of alternating traces. However rather than considering detection of “simple” local maxima/minima as in Definition 4.5, we refer to detection of a generalised notion of local maxima/minima of a trace, that is, maxima and minima which are distanced, at least, by a certain value $\delta$. We formalise this notion in the next definition.

**Definition 4.7** ($\delta$-separated local maxima). Let $\delta \in \mathbb{R}^{+}$, and $\sigma_{A}$ the A projection of a trace $\sigma$ of an $n$-dimensional DESP $\mathcal{D}$ population model. A state $\sigma_{A}[i]$ is the $j$-th, $j \in \mathbb{N}_{>0}$, $\delta$-separated local maximum (minimum), denoted $M_{\delta}[j]$ ($m_{\delta}[i]$), if $i$) it is the largest local maximum (smallest local minimum) whose distance from the preceding $\delta$-separated local minimum $m_{\delta}[j-1]$ (maximum $M_{\delta}[j-1]$) is at least $\delta$ and $ii$) from it, it is possible to reach a smaller (larger) $\delta$-separated state.

To understand the meaning of $\delta$-separated minima and maxima definition let us consider the example in Figure 4.4(b). The nature of the first $\delta$-separated point (either a maximum or a minimum) depends on the initial profile of the trace: if the trace, from its initial state $\sigma_{A}[0]$, first enters the region $\sigma_{A}[0] + \delta$ than the first $\delta$-separated point will be a maximum (conversely if it first enters the $\sigma_{A}[0] - \delta$ region it will be a minimum). For the trace in Figure 4.4(b) the first $\delta$-separated point is the maximum $M_{\delta}[1]$ since none of the preceding local minima $m[1], m[2]$ and local maximum $M[1]$ is sufficiently far apart from the initial state $\sigma_{A}[0]$ (i.e. $m[1], m[2]$ and $M[1]$ have distance less than $\delta$ from $\sigma_{A}[0]$). Observe that the detection of $M_{\delta}[1]$ as the actual first $\delta$-separated maximum is completed only on entering the $(M_{\delta}[1]-\delta)$ region, which happens on observing the downward jump that follows local maximum $M[3]$. To better understand that detection of $\delta$-separated minima/maxima may involve a “temporary detection” phase followed by a “detection finalisation” let us continue the unfolding of the trace. Having entered the $M_{\delta}[1] - \delta$ region we encounter the local minimum $m[4]$, which is temporarily detected as the first $\delta$-separated minimum, since in fact it is distanced more than $\delta$ from $M_{\delta}[1]$. However $m[4]$ is followed by $m[5]$ a smaller local minimum reached without exiting the $(m[4]+\delta)$ region. Thus $m[5]$ replaces $m[4]$ as temporary first $\delta$-separated minimum. $m[5]$ becomes the actual first $\delta$-separated minimum only on entering the $(m[5]+\delta)$ region (which in this case corresponds with reaching of state $M_{\delta}[2]$) since no smaller local minima has been detected before entering $(m[5]+\delta)$.

**Definition 4.8** ($\delta$-separated alternating trace). A trajectory $\sigma$ of an $n$-dimensional DESP $\mathcal{D}$ population model is said $\delta$-separated *alternating* with respect to the $i^{th}$ ($1 \le i \le n$) observed species of $\mathcal{D}$, if it contains infinitely many $\delta$-separated local minima (and equivalently $\delta$-separated local maxima), where $\delta \in \mathbb{R}^{+}$.

For $\sigma_{A}$ a $\delta$-separated alternating trace we denote $\sigma_{A}^{M_{\delta}} = M_{\delta}[1], M_{\delta}[2], \dots$, respectively $\sigma_{A}^{m_{\delta}} = m_{\delta}[1], m_{\delta}[2], \dots$, the projection of $\sigma_{A}$ consisting of the $\delta$-separated local maxima, respectively minima, of $\sigma_{A}$. Observe that the $\delta$-separated max/min (Figure 4.4(b)) are a subset of the “simple” max/min (Figure 4.4(a)). Furthermore the following property holds:
**Proposition 7.** For $\delta=1$ the sequence of $\delta$-separated maxima (minima) of an alternating trace $\sigma_A$ coincides with the list of local maxima (minima), that is: $\sigma_A^{M_1} = \sigma_A^M$ and $\sigma_A^{m_1} = \sigma_A^m$.

The detection of the $\delta$-separated local maxima (minima) for a trace $\sigma_A$ can be described in terms of an iterative procedure through which the list of detected max/min are constructed as $\sigma_A$ unfolds. Such a procedure is formally implemented by the LHA $A_{peaks}$ (Figure 4.5) which we introduce later on. Here, based on the example illustrated in Figure 4.4(b), we informally summarise how detection of $\delta$-separated max/min works. The detection requires storing of the most recent (temporary) $\delta$-separated max (min) into a variable named $x_M(x_m)$, while once detection of a $\delta$-separated maximum (minimum) is completed the corresponding variable $x_M(x_m)$ is copied into a dedicated list, named $Lmax$, resp. $Lmin$ (see Table 4.2), which is filled in with the detected $\delta$-separated max/min points. To understand how detection through automaton $A_{peaks}$ works let us consider again the sample path in Figure 4.4(b). The first element encountered is the local minimum $m[1]$ which is then stored into $x_m = m[1]$. As the trace further unfolds the subsequent maxima (green points) are ignored as long as their distance from the temporary minimum $x_m$ is less than $\delta$, as is the case with $M[1]$. Similarly any local minimum $m[i]$ (yellow point) that is encountered after that stored in $x_m$ is ignored (e.g., $m[2]$), unless it is smaller than $x_m$, in which case $x_m$ is updated with the newly found smaller minimum. As $\sigma_A$ unfolding proceeds we find the next local max $M[2]$ which is distant more than $\delta$ from the temporary minimum $x_m$: this means that $x_m$ currently holds an actual $\delta$-distanced minimum hence its value is appended to $Lmin$ and the procedure starts over, in a symmetric fashion, for the detection of the next maximum.

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

In the remainder we introduce the LHA $A_{peaks}$ which formally implements the detection of the $\delta$-separated peaks of alternating traces.

**The automaton** $A_{peaks}$. We introduce an LHA, denoted $A_{peaks}$ (Figure 4.5), designed for detecting $\delta$-distanced local maxima/minima along alternating traces of a given observed species called $A$. It requires a parameter $\delta$ (the chosen noise level) and the partition of the event set $E = E_{+A} \cup E_{-A} \cup 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
<figure><img src="image_27.png" /><figcaption>Figure 4.5: The $\mathcal{A}_{peaks}$ LHA for detecting local maxima/minima (for observed species A) of noisy periodic traces where local maxima/minima are detected with respect to a chosen level of noise $\delta$.</figcaption></figure>

automaton to loop between **Max** and **Min** hence registering the minima/maxima while doing so. The detailed behavior of $\mathcal{A}_{peaks}$ is as follows. Processing of a trace starts with a configurable filter of the initial transient (represented as a box in Figure 4.5) through which a simulated trace is simply let unfolding for a given $initT$ duration. 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. $x < A - \delta$) or a decrease (i.e. $x > A + \delta$) of the population of the observed species A beyond the chosen level of noise $\delta$. Once within the **Max**→**noisyDec** →**Min**→**noisyInc** loop the detection of local maxima and minima begins. Location **Max** (**Min**) is entered from **noisyInc** (**noisyDec**) each time a sufficiently large (w.r.t. $\delta$) 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 $\delta$ 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]$)$^5$ before storing the new value of $A$ in $x$ and finally increase $n_M$ ($n_m$) the

<footnote>$^5$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</footnote>
<table><caption>Table 4.2: The data variables of automaton $A_{peaks}$ of Figure 4.5 for locating the peaks of a noisy oscillatory traces</caption><thead><tr><th colspan="4">Data variables</th></tr><tr><th>name</th><th>domain</th><th>update definition</th><th>description</th></tr></thead><tbody><tr><td>t</td><td>R<sub>≥0</sub></td><td>reset</td><td>time elapsed since beginning measure (first non-spurious period)</td></tr><tr><td>n<sub>M</sub>(n<sub>m</sub>)</td><td>N</td><td>increment</td><td>counter of detected local maxima (n<sub>M</sub>), minima (n<sub>m</sub>)</td></tr><tr><td>x</td><td>N</td><td>current value of observed species A</td><td>(overloaded) variable storing most recent detected maximum/minumum</td></tr><tr><td>Smax(Smin)</td><td>N</td><td></td><td>sum of detected maxima (minima)</td></tr><tr><td>Lmax<>[Lmin<>]</td><td>N<sup>n</sup></td><td></td><td>array of frequency of heights of detected maxima (minima)</td></tr></tbody></table>

counter of detected maxima (minima). Once in **Max** (**Min**) we stay there as long as we observe the occurrence of reactions which do not decrease (increase) the value of A, hence either reactions of $E_{+A}$ ($E_{-A}$) or of $E_{=A}$. While in **Max** (**Min**) if we observe a reaction of $E_{+A}$ ($E_{-A}$), then we store the new increased (decreased) value of A in x, which becomes the new potential next local maximum (minimum). On the other hand while in **Max** (**Min**) if a “decreasing” (“increasing”) reaction $E_{-A}$ ($E_{+A}$) occurs we move to **noisyDec** (**noisyInc**) from which we can either move back to **Max** (**Min**), if we observe a new increase (decrease) that makes the population of A overpass x (x overpass A), or eventually entering **Min** (**Max**) as soon as the observed decrease (increase) goes beyond the chosen $\delta$ (see above). For the automaton $A_{peaks}$ depicted in Figure 4.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. The data variables of $A_{peaks}$ are summarised in Table 4.2.

**HASL expressions associated with $A_{peaks}$**. We define different HASL expressions to be associated with automaton $A_{peaks}$.

<footnote>given level of the observed population, thus $Lmin_1$ counts the frequency of observed minimum at value 1, $Lmin_2$ the observed minima at value 2 and so on. The number of required $Lmin_i$, $Lmax_j$ variables, which is potentially infinite, can be actually bounded without loss of precision to a sufficiently large value $Lmin_m$ (resp. $Lmax_m$) which must be established manually beforehand, for example by observing few previously generated traces.</footnote>
* $Z_{max} \equiv E[\text{last}(S_{max})/n_M]:$ corresponding to the expected value of the average height of the maximal peaks for the first $N$ detected maxima.

* $Z_{min} \equiv E[\text{last}(S_{min})/n_m]:$ same as $Z_{max}$ but for minima.

* $Z_{PDFmax} \equiv E(\text{last}(L_{max})/n_M):$ enabling to compute the PDF of the height (along a path) of the maximal peaks

* $Z_{PDFmin} \equiv E(\text{last}(L_{min})/n_m):$ enabling to compute the PDF of the height (along a path) of the maximal peaks

Expression $Z_{max}$ ($Z_{min}$) represents the average value of the detected $\delta$-separated maxima (minima). This is obtained by considering the sum of all detected $\delta$-separated local maxima (minima), which is stored in $S_{max}$ ($S_{min}$) and dividing it by the number of detected maxima $n_m$ ($n_m$). Expression $Z_{PDFmax}$ ($Z_{PDFmin}$) allows to estimate the PDF of the height of the detected $\delta$-separated local maxima (minima). This is achieved by dividing the frequency counters of each detected maximal (minimal) peak's height, whose values are stored in array $L_{max}$ ($L_{min}$), by $n_M$ ($n_m$), the number of detected maxima (minima).

**Theorem 3.** If $\sigma_A$ is a $\delta$-separated (with $\delta \in \mathbb{R}^+$) alternating trace then it is accepted by automaton $A_{peaks}$ with parameters $\delta$ and $N \in \mathbb{N}$.

*Proof.* See [Bal15a]. □

### 4.2.4 Case study: circadian clock

In order to demonstrate the HASL based analysis of stochastic oscillator we introduced through the $A_{per}$ et $A_{peaks}$ automata we considered a popular example of oscillator, the so-called *circadian clock*, a biological mechanism responsible for keeping track of daily cycles of light and darkness. Specifically the circadian clock model we considered is that proposed by Vilar *et al.* [VKBL02]: it consists of a genetic circuit responsible for the expression of two co-related proteins schematically depicted in Figure 4.6 and formally expressed by 16 chemical reactions in Equation (4.10). Protein expression is a two steps process: in the first phase a gene transcribes a messenger RNA (mRNA) molecule; in the second phase the mRNA molecule is translated into the target protein. For the model of circadian clock we consider here we denote $M_A$, the mRNA species transcribed by gene $D_A$, and $M_R$ the mRNA transcribed by gene $D_R$. $M_A$ and $M_R$ are then translated into proteins $A$,
<figure><img src="image_28.png" /><figcaption>Figure 4.6: Circadian Clock oscillator network: gene $D_A$ expresses activator protein A through transcription of mRNA $M_A$, while gene $D_R$ expresses the repressor protein R through transcription of mRNA $M_R$.</figcaption></figure>

respectively R.

$$
\begin{align*}
R_1: A + D_A \xrightarrow{\gamma_A} D'_A &\qquad R_9: M_A \xrightarrow{\beta_A} M_A + A \\
R_2: D'_A \xrightarrow{\theta_a} A + D_A &\qquad R_{10}: M_R \xrightarrow{\beta_R} M_R + R \\
R_3: A + D_R \xrightarrow{\gamma_R} D'_R &\qquad R_{11}: A + R \xrightarrow{\gamma_C} C \\
R_4: D'_R \xrightarrow{\theta_R} D_R + A &\qquad R_{12}: C \xrightarrow{\delta_A} R \\
R_5: D'_A \xrightarrow{\alpha'_A} M_A + D'_A &\qquad R_{13}: A \xrightarrow{\delta_A} \emptyset \\
R_6: D_A \xrightarrow{\alpha_A} M_A + D_A &\qquad R_{14}: R \xrightarrow{\delta_R} \emptyset \\
R_7: D'_R \xrightarrow{\alpha'_R} M_R + D'_R &\qquad R_{15}: M_A \xrightarrow{\delta_{M_A}} \emptyset \\
R_8: D_R \xrightarrow{\alpha_R} M_R + D_R &\qquad R_{16}: M_R \xrightarrow{\delta_{M_R}} \emptyset
\end{align*}
\tag{4.7}
$$

Protein A acts as an *activator* for both genes by attaching to promoter region of $D_A$ and $D_R$ (i.e. when A is attached to a gene the mRNA transcription increases). Species $D'_A$ and $D'_R$ represent the state of gene $D_A$, respectively $D_B$, when an activator molecule (A) is attached to their promoter. Note that gene $D_R$ acts as a *repressor* of $D_A$ since when A bounds to its promoter $D_R$ sequesters the activator A and, as a result, the transcription of $D_A$ slows down. The repressing role of $D_R$ is further due to the fact that the expressed protein R inactivates the activator A by binding to it and forming the complex C. Finally the model in Figure 4.6 accounts for degradation of all species: thus the mRNAs $M_A$ and $M_R$, as well as the expressed proteins A and B degrades with given rates (see Table 4.3). Notice that protein A degrades also when attached to R (i.e. when in complex C), and, as a consequence, C turns into R at a rate equivalent to the degradation rate of A.
<table><caption>Table 4.3: reactions' rates for the circadian oscillator</caption><tr><td>α<sub>A</sub></td><td>50 h<sup>-1</sup></td><td>α<sub>R</sub></td><td>0.01 h<sup>-1</sup></td><td>δ<sub>A</sub></td><td>1 h<sup>-1</sup></td><td>δ<sub>R</sub></td><td>0.2 h<sup>-1</sup></td></tr><tr><td>α'<sub>A</sub></td><td>500 h<sup>-1</sup></td><td>α'<sub>R</sub></td><td>50 h<sup>-1</sup></td><td>γ<sub>A</sub> = γ<sub>R</sub></td><td>1 mol<sup>-1</sup>h<sup>-1</sup></td><td>γ<sub>C</sub></td><td>2 mol<sup>-1</sup>h<sup>-1</sup></td></tr><tr><td>β<sub>A</sub></td><td>50 h<sup>-1</sup></td><td>β<sub>R</sub></td><td>5 h<sup>-1</sup></td><td>θ<sub>A</sub></td><td>50 h<sup>-1</sup></td><td>θ<sub>R</sub></td><td>100 h<sup>-1</sup></td></tr><tr><td>δ<sub>MA</sub></td><td>10 h<sup>-1</sup></td><td>δ<sub>MR</sub></td><td>0.5 h<sup>-1</sup></td><td></td><td></td><td></td><td></td></tr></table>

<figure><img src="image_29.png" /><figcaption>Figure 4.7: GSPN encoding of the system (4.10) of chemical equations corresponding to the circadian-clock.</figcaption></figure>

The model of Figure 4.6 corresponds to the system of chemical equations (4.10), whose (continuous) kinetic rates (taken from [VKBL02]) are given in Table 4.3.

**Stochastic model** Equations (4.10) can give rise to either a system of ODEs or to a stochastic process. Here we focus on the discrete-stochastic semantics: Figure 4.7 shows the GSPN encoding of equations (4.10) developed with COSMOS. The configuration of the GSPN (i.e. the stochastic process) requires setting the initial population and the rates of each transition (i.e. reaction). For the initial population, following [VKBL02], we observe that the model comprises one gene $D_A$ and one $D_R$, which can either be in free-state (no activator $A$ is attached to the promoter) or in activator-bound state, i.e. $D'_A$, respectively $D'_R$. As a consequence the population of species $D_A$ and $D_R$ is bounded by the following invariant constraints: $D_A + D'_A = 1$ and $D_R + D'_R = 1$ (in fact places $DA$, $DA'$ and $DR$, $DR'$ of net in Figure 4.7 are the only places covered by P-invariants). The remaining species are initially supposed to be “empty”, hence they are initialised to 0. Concerning the transition rates, for simplicity we assume a unitary volume of the system under consideration, hence all continuous rates in Table 4.3 can be used straightforwardly as rates of the corresponding discrete-stochastic reactions. In this case we assume all reactions following a negative exponential law.

The oscillatory dynamics of the GSPN model of Figure 4.7 can be observed by plotting of a simulated trajectory (Figure 4.8). Observe that the frequency of oscillations varies considerably with the degradation rate of the repressor ($R$) protein: a faster degradation
<figure><img src="image_30.png" /><figcaption>Figure 4.8: Single trajectory showing the oscillatory character of activator *A* and repressor *B* dynamics with normal repressor’s degradation rate $\delta_R$=0.2 (left) and with 10× speed-up, i.e. $\delta_R$ = 2 (right).</figcaption></figure>

of *R* (right), intuitively, results in a higher frequency of oscillations. In the remainder we formally assess the oscillatory characteristics (i.e. the period and the peaks of oscillations) of the circadian clock model by application of the previously described approach, i.e. by analysing the stochastic process deriving from synchronisation of the circadian clock GSPN model with the *A*<sub>per</sub> and *A*<sub>peaks</sub> automata.

**Measuring the period of the circadian clock.** We performed a number of experiments aimed at assessing the effect that the degradation rate of the repressor protein ($\delta_R$) has on the period of the circadian oscillator. Figure 4.9 (right) shows three plots representing the PDF of the period (obtained through the HASL formula ($A_{per}, PDF(Last(t)/N$)) for three values of $\delta_R$. With $\delta_R$ = 0.2 (i.e. the original value as given in [VKBL02]) the PDF is centred at $t$ = 24.9, i.e. slightly more of the standard 24 hours period expected for a circadian clock. On the other hand speeding up the repressor degradation of 10 times (i.e. $\delta_R$ = 2) yields a slightly more than halved oscillation period (i.e. PDF centred at $T$ = 10.8). Finally slowing down the degradation rate of a half (i.e. $\delta_R$ = 0.1) yields a less than doubled oscillation period (i.e. PDF centred at $T$ = 40.7). Notice that plots in Figure 4.9 (right) also indicate that slowing down of the degradation of the repressor *R* affects as well the variability of the oscillation period as witnessed by the increasing width of the bell-shape form of the PDF plots in Figure 4.9 (right).

Figure 4.9 (left) shows plots for the period mean value (red plot) and the period variance (blue plot) obtained through by iterative evaluation of the HASL formulae ($A_{per}, E[Last(\bar{t}_p)]$), respectively, ($A_{per}, E[Last(s_{tp}^2)]$) in function of the degradation rate $\delta_R$. They indicate that slowing down the degradation of the repressor yields, on one hand, to a lower the period of oscillations, and on the other, augmenting the irregularity of the periods (i.e. augmenting the period’s variances). All plots in Figure 4.9 result from sampling of
<figure><img src="image_31.png" /><figcaption>Figure 4.9: The PDF (left) and the mean value vs the variance (right) of the period of oscillations of protein *A* of the circadian clock measured with *A*<sub>per</sub> in function of the repressor’s degradation rate.</figcaption></figure>

finite trajectories consisting of *N* = 100 periods, where periods have been detected using *L* = 1 and *H* = 1000 as partition thresholds, and target estimates have been computed with confidence level 99 and confidence-interval width of 0.01. Furthermore the PDF plots in Figure 4.9 (right) have been computed using a discretisation of the period support interval [0, 50] into subintervals of width 0.1.

**Measuring the peaks of oscillations of the circadian clock.** We performed a num-
ber of experiments aimed at assessing the effect that the degradation rate of the repressor
protein ($\delta_R$) has on the peaks of oscillation for both protein *A* and *R*. Figure 4.10 shows
plots for the mean value of the minimal and maximal peaks of oscillations for both *A* and
*R*. Results indicate that while the degradation rate $\delta_R$ has no effect on the oscillation
peaks of *A* (both maximal and minimal peaks of *A* are constant independently of $\delta_R$), it
affects the maximal peaks (only) of *R*. Specifically the mean value of *R*'s maximal peaks
decreases with the increasing of $\delta_R$ (while the minimal peaks of *R* are constantly at 0),
notice that this is in agreement with what indicated by the single trajectories depicted
in Figure 4.8, as these show that the oscillations of both *A* and *R* consist of positive
half-period (where the population is positive) interleaved by zero half-period (where the
population is zero). Notice that Figure 4.10 contains also plot for the absolute maximum
of population of *A* and *R* measured along the sampled trajectories through trivial HASL
expressions $Z \equiv AVG(max(x))$ (where *x* is a variable used to record the population of the
observed species along a synchronising path). All plots in Figure 4.10 result from sampling
of finite trajectories containing of *N* = 100 maximal peaks and using a noise parameter
$\delta$ = 10% AVG(max(x)), meaning that for evaluating the mean value of maximal peaks we
discarded all critical points (i.e. the noise spikes) distanced one another less than 10%
<figure><img src="image_32.png" /><figcaption>Figure 4.10: The mean value of the minimal and maximal peaks of proteins A and R of the circadian clock measured with $A_{peaks}$ in function of the repressor’s degradation rate.</figcaption></figure>

of the absolute maximum of the observed species. Finally, again points of every plot in Figure 4.10 have been computed with confidence level 99 and confidence-interval width of 0.01.

## 4.3 Tuning stochastic oscillators

The second contribution described in this chapter is that of a framework for tuning of stochastic oscillators, that is, a framework that allows for searching the parameter space of an oscillator model so to identify the parameter’s values that comply with a target oscillation period. Such framework combines the partition based HASL characterisation of noisy periodicity outlined in the first part of this chapter with the Automaton-ABC parameter inference framework outlined in Chapter 3. As we have seen the adaptation of the ABC parameter inference scheme to HASL relies on the existence of a notion of distance between the paths sampled from a given instance of a model and the sought (target) behaviour. Therefore, in this case, we need to identify a suitable distance that can effectively drive the inference of parameters according to a desired oscillation character, in this case the oscillation period.

**Definition 4.9** (distance from target period). For $\sigma_A$ a noisy periodic trajectory and $\bar{t}_p^{(obs)} \in \mathbb{R}_{>0}$ a target mean period duration we define the distance of $\sigma_A$ from $\bar{t}_p^{(obs)}$ w.r.t. the first $n \in \mathbb{N}$ period realisations as

$$ \mathrm{dist}(\sigma_A, n, \bar{t}_p^{(obs)}) = \mathrm{dist}(\bar{t}_p(n), s_{\bar{t}_p}^2(n), \bar{t}_p^{(obs)}) = \min(\frac{\bar{t}_p(n) - \bar{t}_p^{(obs)}}{\bar{t}_p^{(obs)}}, \frac{\sqrt{s_{\bar{t}_p}^2(n)}}{\bar{t}_p^{(obs)}}) \quad (4.8) $$
where $\bar{t}_p(n)$ ($s_{\bar{t}_p}^2(n)$) denotes the sample mean value (the variance) computed w.r.t. the first $n$ periods detected along $\sigma_A$ (as per Definition 4.2 and Definition 4.4).

Notice that distance (4.8) establishes a form of multi-criteria selection of parameters as both the mean value and the variance of the detected periods are constrained. For example with a 10% tolerance (i.e. $\epsilon = 0.1$ in ABC terms) only the parameters $\theta \in \Theta$ that issue a relative error (w.r.t. the target mean period) not beyond 0.1 are selected.

### 4.3.1 An automaton for the distance from a target period

In Figure 4.11 we introduce the automaton $\mathcal{A}_{per}^{\bar{t}_p^{(obs)}}$, i.e, an adaptation of the $\mathcal{A}_{per}$ (Figure 4.3), suitable for assessing the distance (as per Definition 4.9) between the mean period measured on paths of DESP oscillator $M_\theta$ and a target period duration $\bar{t}_p^{(obs)}$. The au-

<figure><img src="image_33.png" /><figcaption>Figure 4.11: $\mathcal{A}_{per}$: an LHA for selecting noisy periodic traces (with respect to an observed species $A$) related to partition $low = (-\infty, L]$, $mid = (L, H)$ and $high = [H, +\infty)$.</figcaption></figure>

tomaton $\mathcal{A}_{per}^{\bar{t}_p^{(obs)}}$ differs from $\mathcal{A}_{per}$ only in the autonomous arcs that leads to the accepting location **end**: differently from $\mathcal{A}_{per}$ those arcs are equipped with an update through which a dedicated extra variable, named $d$, is set with the value of the measured distance as per Equation (4.8) (see Table 4.4 for the list of variables).

### 4.3.2 Case studies

We demonstrate the HASL-ABC for tuning stochastic oscillator through a couple of case studies.
<table><caption>Table 4.4: Variables of the $A_{per}^{\bar{t}_p^{(obs)}}$ automaton.</caption><thead><tr><th>name</th><th>domain</th><th>update definition</th><th>description</th></tr></thead><tbody><tr><td>n</td><td>N</td><td>increment</td><td>counter of detected periods</td></tr><tr><td>t<sub>p</sub></td><td>R≥0</td><td>reset</td><td>duration last period</td></tr><tr><td>t&#773;<sub>p</sub></td><td>R≥0</td><td>f(t&#773;<sub>p</sub>, t<sub>p</sub>, n) = (t<sub>p,n</sub>·n+t<sub>p</sub>) / (n+1)</td><td>period mean</td></tr><tr><td>s<sup>2</sup><sub>t<sub>p</sub></sub></td><td>R≥0</td><td>g(s<sup>2</sup><sub>t<sub>p</sub></sub>, t&#773;<sub>p</sub>, t<sub>p</sub>, n) = (n-1)/(n-2) · s<sup>2</sup><sub>t<sub>p</sub></sub> + ((t<sub>p</sub>-t<sub>p</sub>)<sup>2</sup>) / n</td><td>period variance</td></tr><tr><td>d</td><td>R≥0</td><td>min(t&#773;<sub>p</sub>-t<sub>p</sub>(obs) / t<sub>p</sub>(obs), √s<sup>2</sup><sub>t<sub>p</sub></sub> / t<sub>p</sub>(obs))</td><td>distance from target period</td></tr></tbody></table>

## A synthetic 3-ways oscillator

The CRN in (4.9) represents a simple model of sustained oscillator called doping 3-way oscillator [Car09].

$$ 
\begin{array}{l}
R_1: A+B \xrightarrow{r_A} 2B \qquad R_2: B+C \xrightarrow{r_B} 2C \qquad R_3: C+A \xrightarrow{r_C} 2A \\
R_4: D_A+C \xrightarrow{r_C} D_A+A \qquad R_5: D_B+A \xrightarrow{r_A} D_B+B \qquad R_6: D_C+B \xrightarrow{r_B} D_C+C
\end{array}
\qquad (4.9) $$

<figure><img src="image_34.png" /></figure>

(a) Synchronisation of a path of the 3-way oscillator with $A_{per}$: coloured segments represent the corresponding $A_{per}$ location, while $n$ is the period counter variable.

<figure><img src="image_35.png" /></figure>

(b) Posterior distribution for the $r_A$ parameter obtained by application of the automaton-ABC with $A_{per}$.

Figure 4.12: A sampled trajectory of the 3-way oscillatory (left) with coloured highlights of $A_{per}$ corresponding location transitions and the resulting posterior distribution (right) obtained for $r_A$ while keeping $r_B=r_C=1$.

It consists of 3 main species A, B and C forming a positive feedback loop (through reactions $R_1, R_2, R_3$) plus 3 corresponding invariant (doping) species $D_A, D_B, D_C$ whose goal is to avoid extinction of the main species (through reactions $R_4, R_5, R_6$). It can be shown that the total population is invariant and that the model yields sustained noisy oscillation for the 3 main species, whose period and amplitude depend on the model's
parameters $r_A$, $r_B$ and $r_C$ (as well as on the total population). Figure 4.12(a) depicts a path sampled from (a given configuration of) model (4.9).

**Experiment 1 (one-dimensional).** In this experiment we considered $s_0 = (A_0, B_0, C_0, (D_A)_0, (D_B)_0, (D_C)_0) = (333, 333, 333, 10, 10, 10)$ as initial state, we fixed the rate constants $r_B = r_C = 1.0$, and estimated the posterior distribution for $r_A$ considering a uniform $U(0, 10)$ prior and a target mean period $\bar{t}_p^{(\text{obs})} = 0.01$. For the automaton $A_{\text{per}}$ the noisy-period dependent partition we considered is $L = 300$ and $H = 360$ while for each trajectory we observed $N = 4$ periods. For the ABC algorithm we used $N = 1000$ particles and considered a 20% tolerance ($\epsilon = 0.2$). Figure 4.12(b) shows the resulting automaton-ABC posterior. We observe that 1) the posterior support being included in the prior is $[0.0, 4.0] \subset [0, 10.0]$, i.e., we have reduced the parameter space to a subset where it is probable to obtain trajectories with a mean period of 0.01 (relative to a 20% tolerance) and 2) the posterior has only one mode, which is quite sensible, as having fixed $r_B$ and $r_C$, one would expect the mean period duration being directly linked to the kinetics of reaction $R_1$, hence to its only parameter, the kinetic rate constant $r_A$.

<figure><img src="image_36.png" /><figcaption>Figure 4.13: Correlation plot of automaton-ABC posterior with $A_{\text{per}}$ posterior for the 3D experiment of doping 3-way oscillator.</figcaption></figure>

**Experiment 2 (three-dimensional).** Here we performed a 3-dimensional version of the previous experiment assuming the same Uniform prior distribution for all 3 parameters, that is, $r_A, r_B, r_C \sim U(0, 10)$, however considering a different target mean period $\bar{t}_p^{(\text{obs})} = 20$. Figure 4.13 shows the correlation plot of the resulting posterior distribution obtained by application of the automaton-ABC method with $A_{\text{per}}$. The one-dimensional histogram on the diagonal of Figure 4.13 represents the marginal distribution of each parameter while the plots in the upper triangular submatrix show the histograms of the two-dimensional marginal distribution (e.g., plot on position (1,3) refers to the marginal
distribution $p(r_A, r_C | \bar{t}_p^{(obs)})$), whereas the lower triangular matrix show the scatter plots of the two-dimensional marginal distributions. Based on the diagonal plots we observe that most of the parameters that results in a period close to the target one ($\bar{t}_p^{(obs)} = 20$) are within the support $[0, 4] \times [0, 3] \times [0, 4]$ and, furthermore the particles form a 3D parabolic shape since the three two-dimensional projections have a parabolic shape, according to plots in the upper triangular part of Figure 4.13. Also, one can notice that for each two-dimensional histogram, the region near point (1, 1) is a high probability one, which is consistent with the previous one-dimensional experiment.

## Repressilator

We consider a model (4.10) of a synthetic genetic network known as Repressilator which was developed to reproduce oscillatory behaviours within a cell [EL00]. It consists of 3 proteins $P_1, P_2, P_3$ forming a negative feedback loop, with $P_1$ repressing $P_2$'s transcription gene $G_2, P_2$ repressing $P_3$'s transcription and so on.

$$ 
\begin{array}{l@{\quad}l@{\quad}l}
\text{R1: } G_1 \xrightarrow{r_1} G_1 + M_1 & \text{R2: } G_2 \xrightarrow{r_2} G_2 + M_2 & \text{R3: } G_3 \xrightarrow{r_3} G_3 + M_3 \\
\text{R4: } M_1 \xrightarrow{\beta} M_1 + P_1 & \text{R5: } M_2 \xrightarrow{\beta} M_2 + P_2 & \text{R6: } M_3 \xrightarrow{\beta} M_3 + P_3 \qquad (4.10) \\
\text{R7: } M_1 \xrightarrow{1} \emptyset & \text{R8: } M_2 \xrightarrow{1} \emptyset & \text{R9: } M_3 \xrightarrow{1} \emptyset \\
\text{R10: } P_1 \xrightarrow{1} \emptyset & \text{R11: } P_2 \xrightarrow{1} \emptyset & \text{R12: } P_3 \xrightarrow{1} \emptyset
\end{array}
$$

Following [EL00] we assumed mass-action dynamics with a common rate constant $\beta$ for translations reactions ($R_4, R_5, R_6$) and constant rate 1 for species degradation ($R_7$ to $R_{12}$). Conversely transcription dynamics ($R_1, R_2, R_3$) is assumed to follow a Hill function dynamics given by the follow parameter definitions $r_1 = \frac{\alpha}{1+\frac{\alpha}{[P_3]^n}} + \alpha_0$, $r_2 = \frac{\alpha}{1+\frac{\alpha}{[P_1]^n}} + \alpha_0$ and $r_3 = \frac{\alpha}{1+\frac{\alpha}{[P_2]^n}} + \alpha_0$, where $n$ is the Hill coefficient, $\alpha$ is related to transcription growth and $\alpha_0$ is the parameter related to the minimum level of transcription growth. The parameter space is therefore 4-dimensional with $\theta = (\alpha, \beta, n, \alpha_0) \in \mathbb{R}^4$.

<figure><img src="image_37.png" /><figcaption>Figure 4.14: $P_1$-projections of paths simulated with $(\beta, n, \alpha_0) = (2, 2, 0)$ and $\alpha \in \{50, 200, 1000\}$.</figcaption></figure>

Each parameter affects the resulting oscillation, e.g., Figure 4.14 shows examples of noisy periodic paths corresponding to different values of $\alpha$, with $(\beta, n, \alpha_0) = (2, 2, 0)$.
<figure><img src="image_38.png" /><figcaption>Figure 4.15: Correlation plot of automaton-ABC posterior with *Aper* posterior for the 3D experiment of repressilator model.</figcaption></figure>

**Experiment 1.** This is a 3-dimensional experiment in which we considered $s_0 = ((M_1)_0, (M_2)_0, (M_3)_0, (P_1)_0, (P_2)_0, (P_3)_0) = (0,0,0,5,0,15)$ as initial state, we fixed $\alpha_0 = 0$, and estimated the posterior distribution for the remaining parameters considering the following prior $\alpha \sim U(50, 5000)$, $\beta \sim U(0.1, 5.0)$, $n \sim U(0.5, 5.0)$. The target mean period was $\bar{t}_p^{(obs)} = 20$, $L = 50$ and $H = 200$ the noisy-period dependent partition. For the ABC algorithm we used $N = 1000$ particles and considered a 10% tolerance ($\epsilon = 0.1$). Figure 4.15 shows the correlation plot of the resulting automaton-ABC posterior. We observe that with this setting the Repressilator oscillations are most sensitive to parameter $n$ whose marginal posterior is much narrower than that of $\alpha$ and $\beta$ (i.e. varying $n$ induces more instability than $\alpha$ and $\beta$).

**Experiment 2.** This is a 4-dimensional version of the previous experiment in which we considered also a uniform prior for $\alpha_0 \sim U(0.1, 5.0)$. Figure 4.16 shows the correlation plot of the resulting automaton-ABC posterior. We observe that adding a degree of freedom on $\alpha_0$ has changed the correlation between $\alpha$ and $\beta$ as well as $\alpha$ and $n$ whereas correlation between $\beta$ and $n$ seems to have the same shape.

## 4.4 Perspectives

As a first contribution in this chapter (Section 4.2) we tackled the problem of using an automata-based formalism to formally analyse stochastic oscillators. In this respect the contribution is twofold. Firstly we have generalised Spieler's time-automata noisy periodicity detection, and shown that through HASL one can assess oscillation related indicators that could not be accounted for otherwise. Secondly we have outlined an alternative
<figure><img src="image_39.png" /><figcaption>Figure 4.16: Correlation plot of automaton-ABC posterior with *Aper* posterior for the 4D experiment of repressilator model.</figcaption></figure>

automata-based approach for analysing periodicity which is aimed at identifying the peaks of stochastic oscillations leading to the characterisation of amplitude related indicators. Although the proposed HASL-based framework opens up to gaining detailed insights about the characteristics of a stochastic oscillators, it provides little help for establishing whether the considered system oscillate sustainably or not. As shown by Spieler the existence of sustained oscillation can be established by checking of a steady-state property on the product Markov chain$^6$, i.e. the ergodic-by-construction Markov chain given by the product of the considered model and the DTA period detector. An interesting perspective would be to investigate whether the HASL framework for oscillation analysis could be adapted to establishing whether the oscillations are sustained. Intuitively that would entail the integration, within the HASL model checking framework, of simulation based methods for the verification of steady-state measures, such as *perfect sampling* [RP09], which indeed is an aspect I am very much willing to investigate. Another interesting evolution of this contribution would be to move towards a full automatisation of periodicity detection. In fact the partition based period detection approach depend upon the (manually) chosen constant parameters *L* and *H* that characterise the partition and therefore both the period detection and the measured indicators (period mean value and variance) are very much

<footnote>$^6$the model oscillate sustainably if converging paths of the product process have zero probability at steady-state</footnote>
affected by the chosen *L* and *H* values. As a consequence if one chooses too "extreme" *L* and *H* values the detection of a period may become a very unlikely (possibly a zero probability) event, therefore, in practice, one must have an idea about the amplitude (and location) of the high and low peaks of the oscillation, in order to correctly set the *L* and *H* parameter of the detector automata beforehand. A partition independent detection of the oscillation periods, intuitively, could be worked out following the peaks based period detection approach outlined in Section 4.2.3, although the integration of partition-based with peak-based period detection has not yet been investigated so far.

As a second contribution of this chapter (Section 4.3) we have shown how to adapt
the ABC-based parametric verification framework outlined in Chapter 3 to the problem
of tuning of stochastic oscillators, that is, to identifying the probability distribution with
which the parameters of a stochastic oscillator match a given period of the oscillations.
Such contribution currently only take into account the oscillation period as the constraining
factor driving the parameter search. An interesting perspective would be to extend this
work so that also the oscillation amplitude can be accounted for during the parameter search
process, something that intuitively could be achieved based on the amplitude related HASL
specifications outlined in Section 4.2.
# Chapter 5

## Data-driven predictive modeling of periodic phenomena

As argued in the previous chapter many real life phenomena exhibit some form of periodicity, examples of which can be found in different domains e.g., road traffic, patient arrival to an emergency department or photovoltaic energy production. In Chapter 4 we tackled the problem of how to automatically assess (probabilistic) periodicity related indicators of a formally given stochastic oscillator model and also how to tune the parameters of such a model so that it exhibits a desired periodicity. In this chapter we consider a complementary aspect, namely, the derivation of (stochastic) models based on observations of periodic phenomena. We stress that the ability to derive models that accurately reproduce the periodicity exhibited by data samples issued by periodic systems is vital for predictions, planning and finally optimisation of such systems. Here we report on a contribution in the context of optimisation of hospitals’ emergency department (ED), for which periodicity is indeed relevant as evidenced by datasets we had access to.

**Outline of the contribution.** In collaboration with Andras Horváth, Roberto Aringhieri and in the context of the PhD thesis work of Davide Duma (all three of them from University of Turin) we studied the problem of deriving a stochastic model of the patient arrivals (to an ED) that is capable of correctly reproducing the statistical character of a given arrivals dataset. To this aim we considered a real dataset relative to the activity of the ED of an hospital of the city of Cantù in northern Italy and which contains different data including the complete set of timestamps of the arrival of each patients over one year (i.e. the year 2015). Based on the statistical analysis of the dataset, which evidenced the periodic character of the arrival process, we derived different instances of Markovian models belonging to different classes (i.e. Markov renewal processes, Markov arrival processes, hidden Markov models) and showed that only one family of model, amongst those considered, is capable of fully reproducing the periodic statistical character of the dataset. To
further validate the derived models we developed a formal encoding in terms of stochastic Petri nets. This allowed us to compare the performances of each arrival models, coupled with a model of ED services, through assessment of a key performance indicator (KPI) formally encoded and accurately assessed through HASL statistical model checking. The obtained results evidenced that one amongst the derived arrival models faithfully match the dataset behaviour hence it can reliably be used for analysing the performances of different ED services.

**Structure of the chapter.** The chapter is organised as follows: Section 5.1 reports on the results of statistical analysis of the ED dataset, evidencing the periodic nature of patients arrivals. In Section 5.2 we provide a comparative analysis of how different families of stochastic models, i.e., Markov renewal processes, Markov arrival processes and hidden Markov model, derived from dataset are capable of adequately reproducing the periodic nature of the observations. Finally in Chapter 5.3 we tackled the validation of derived models. To this aim we considered a simple model of ED services (representing the services a ED patient undergoes after triage), we combined it with each derived model as well as a with a model reproducing exactly the dataset arrivals and by assessing relevant performance indicators on the arrivals/service combined models we established that only one kind of derived model is capable of adequately matching the dataset.

## 5.1 Statistical analysis of an emergency department dataset

The first goal of our study was to gain an understanding on relevant characteristics of the patient's arrival process by analysing the one year timespan dataset we had access to. We considered two perspectives: that of the patients' *inter-arrival time* and that of the patients' *number of arrivals per time frame*, which we represented through different families of random variables given in the following definition.

**Definition 5.1.** For $D$ a dataset spanning over several weeks time frame and whose entries are equipped with timestamps relative to the occurrence of a given type of event we consider the following families of random variables:

* $X_i$ ($i = 1, 2, 3, ...$): the $i$-th inter-arrival time, i.e., the time elapsed between the $i$-th and $(i+1)$-th occurrence of the considered kind of event.

* $Y_{H,i}$, $Y_{D,i}$ and $Y_{W,i}$ ($i = 1, 2, 3, ...$): the number of occurrences of the considered event during the $i$-th hour ($Y_{H,i}$), day ($Y_{D,i}$), respectively week ($Y_{W,i}$).

Clearly, variables of families $Y_{H,i}$, $Y_{D,i}$ and $Y_{W,i}$ are less informative than those of family $X_i$ and yet they are relevant as they allow us for easily investigate the periodic nature of the ED's patient arrival process. For example, the variables $Y_{H,3+24j}$ with $j = 0, 1, 2, ...$
<figure><img src="image_40.png" /><figcaption>Figure 5.1: Mean and variance of number of arrivals in the i-th hour of the day.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 5.2: Pmf of number of arrivals between 1pm and 2pm (left) and between 5am and 6am (right) and pmf of the Poisson distribution with equal mean.</figcaption></figure>

<figure><img src="image_42.png" /></figure>

refer to the number of arrivals between 2am and 3am during the 1st, the 2nd, the 3rd day, and so on.

For both families of variables we computed statistical indicators such as the sample mean, the sample variance, the sample probability distribution and the sample autocorrelation function, the latter aimed at investigating the periodic character of the patient arrival process.

**Arrival per hour of the day.** Based on variables of the $Y_{H,i}$ family we estimated the mean value $E_{H,i} = E[Y_{H,i+24j}|j=0,1,2,...]$ and the variance $Var_{H,i} = Var[Y_{H,i+24j}|j=0,1,2,...]$ of the number of patient arrivals during the i-th hour of the day which are both depicted in Figure 5.1 in function of i. The highest mean number of patient arrivals is about 5.9 per hour and is registered between 9am and 10am whereas the lowest number of arrivals is instead less than 1 patient per hour and takes place between 5am and 6am: as expected, the time of the day has a strong impact on the number of patient arrivals. We observe that the variance, $Var_{H,i}$, has very similar values to those of the mean $E_{H,i}$, which is compatible with a well-known characteristic of the Poisson process, hence suggesting that the patient arrival process may be adequately modelled through a Poisson process with time inhomogeneous intensity. To further investigate this aspect, we compare in Figure 5.2 the probability distribution of the number of arrivals during a given hour of the day computed over the dataset with the corresponding Poisson distribution i.e. that with the same mean<fnref n="1" />. Plots in Figure 5.2 compare the resulting probability mass function (pmf) for two different time intervals: 1pm-2pm which is the hour of the day where the Poisson distribution differs most<fnref n="2" /> from the actual distribution of the number of arrivals and the interval 5am-6am is the hour where the Poisson distribution is most similar to the experimental distribution.

To investigate the periodicity of one-hour intervals we considered the hourly autocor-

<footnote n="1">¹the Poisson distribution with mean given by the $E_{H,i}$ for the corresponding hour i.</footnote>

<footnote n="2">²the difference is calculated as the sum of absolute differences in the pmf.</footnote>
<figure><img src="image_43.png" /><figcaption>Figure 5.3: $R_{H,n}$, i.e., acf of the sequence $Y_{H,i}$ as function of $n$.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 5.4: Average and variance of number of arrivals in the ith day of the week.</figcaption></figure>

relation function (acf) defined as:

$$R_{H,n} = \frac{E[(Y_{H,i} - E[Y_{H,i}]) (Y_{H,i+n} - E[Y_{H,i}])]}{Var[Y_{H,i}]} \quad \text{with } n = 1, 2, 3, \dots$$

and whose value is depicted in Figure 5.3 in function of $n$. The persistent periodicity of the hourly arrivals is proved by the sustained oscillation exhibited by the autocorrelation plot. The periodic nature of the autocorrelation plot is justified by the non-constant hourly patient arrival intensity (as proved by the mean value plot in Figure 5.1) and its sustained character is evidenced by the fact that autocorrelation oscillations remain unaltered even for very large values of $n$.

**Arrival per day of the week.** To figure out whether the periodicity is present also w.r.t. the day of the week we performed the same type of analysis but taking into account variables $Y_{D,i}$. In Figure 5.4 we depict $E_{D,i} = E[Y_{D,i+7j}|j = 0, 1, 2, \dots]$ and $Var_{D,i} = Var[Y_{D,i+7j}|j = 0, 1, 2, \dots]$ with $i = 1, 2, \dots, 7$ and observe that the mean number of arrivals is essentially constant over the days of the week whereas the variance differs to a large extent from the mean on most days of the week (which makes the daily arrival process incompatible with a Poisson process). Furthermore the autocorrelation calculated also w.r.t. $Y_{D,i}$ (not shown here) exhibited only a very mild correlation between number of arrivals of consecutive days therefore we may assert the absence of periodicity when considering arrivals per day of the week, which indicates we will not have to account for it when developing models aimed at reproducing the statistical character of the dataset.

**Inter-arrival analysis.** Finally, we considered the family of random variables $X_i$ ($i = 1, 2, 3, \dots$) representing the patients’ inter-arrival times and computed the same statistical indicators. The resulting mean inter-arrival time is $E[X_i] = 19.43$ (minutes), the variance $Var[X_i] = 483.18$, while the shape of the probability density function computed over $X_i$ is depicted in Figure 5.5. Figure 5.6 instead depicts the autocorrelation computed for $X_i$. The damped oscillation profile exhibited by the autocorrelation is due to the fact that
<figure><img src="image_45.png" /><figcaption>Figure 5.5: Front (left) and tail (right) of the pdf of the inter-arrival times.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Figure 5.6: Autocorrelation function of the inter-arrival times.</figcaption></figure>

although periodicity exists between the inter-arrivals The autocorrelation present in $X_i$ is depicted in Figure 5.6 where there is oscillation due to the fact that a period with high arrival intensity (8am-22pm) are followed by a period with low arrival intensity (22pm-8am). This oscillation however is less easy to interpret than that in Figure 5.3 and fades away as $n$ grows because for large values of $n$ there is no deterministic connection between the time of the day of the $i$th and the $(i + n)$th arrival.

## 5.2 Deriving adequate models of patient arrivals

We seek to derive adequate stochastic models capable of reproducing the statistical character of a dataset (in our case the Cantu hospital ED patients' arrival dataset) including the sustained periodicity of the number of arrivals per hour it exhibits. The search for an adequate model involves selecting a suitable family of stochastic processes and, by application of a parameter estimation approach, identifying a model instance that adequately approximate the dataset. The choice of the target family of stochastic process depends very much on the statistical perspective one wishes to mimic.

For example to adequately reproduce the PDF of the inter-arrival time we may consider the family of Markov renewal processes (MRPs) and derive the parameters of the corresponding continuous phase-type (CPH) distribution. Although MRPs can be used to obtain an arbitrarily precise approximation of the inter arrival time PDF they fall short when it comes to reproduce periodicity as they do not allow to account for existing correlation in the observed inter arrival times, as is indeed the case with the considered dataset (Figure 5.6). In this respect the family of Markov arrival processes (MAPs) constitute an improvement as they do allow one to account for correlation between inter-arrival times hence they allow for generating arrivals that not only match the inter-arrival PDF but that are also correlated. However if the considered dataset exhibits correlation also w.r.t. the number of arrivals per time-frame (as in our case), MAPs may not be a convenient option to mimic both kind of correlation. In fact the ability of MAPs to cope with correlated
data affects the size of the the model (i.e. the number of phases it consists of). In the
remainder we show that the size of a MAP capable of reproducing (to some extent) the
correlation between inter-arrivals times is already critical, hence if we wanted to reproduce
also the correlation between the number of arrivals per hour we would need an even larger
MAP for which, however, fitting techniques would not be applicable due to the too large
number of states. Therefore in order to obtain a model capable of matching the hourly
periodicity of the arrival process, we resort to the class of hidden Markov models (HMMs)
and we apply *learning approaches* to derive models' instances from the dataset.

5.2.1 Markovian renewal process approximation

The first family of models we consider for approximating the arrival process is that of Markov renewal process, i.e. a kind of process suitable for modelling independent and identically distributed inter-arrival times. With a MRP the dataset inter-arrival time distribution is approximated by a degree-n phase type (PH) distribution, i.e. the distribution of the time to absorption of a (continuous-time) Markov chain consisting of n transient states (the phases) and one absorbing state [Neu75]. A degree-n CPH distribution is conveniently described by the vector $\alpha \in [0, 1]^n$ of initial probabilities of the transient states<fnref n="3" /> and the matrix $Q \in \mathbb{R}^{n \times n}$ consisting of the transition intensities among the transient states and the opposite of the sum of the intensities of the outgoing transitions in the diagonal (this allows to determine the transition intensities toward the absorbing state). Based on basic properties of CTMCs the analytical expression of the pdf of the time to absorption $f: \mathbb{R}^+ \rightarrow [0, 1]$ of a CPH is given by $f(x) = \alpha e^{Qx}(-Q)I$ where $I$ is the column vector of 1s, and, consequently, that of its i-th ($i \in \mathbb{N}^+$) moment is given by $m_i = i!\alpha(-Q)^{-i}I$.

Figure 5.7 illustrates a toy example of a degree-3 CPH whose initial probability vector $\alpha$ and intensities matrix $A$ are given in the left, the corresponding state-transition graph is depicted in the middle (the probability values inside the states indicates the initial probability of the state while real values labelling the arcs represent transition intensities and the gray state is the absorbing one) and the corresponding PDF of the time to reach the absorbing state is shown on the right.

Having fixed the degree $n$ for the target CPH, our goal is to derive $\alpha$ and $Q$ such that the pdf of the associated CPH distribution is a good approximation of the patient inter-arrivals pdf. There are two families of approaches to face this problem, i.e., those based on the maximum likelihood method as opposed to those based on moments matching methods. Maximum likelihood methods are computationally more hungry as they require to consider the entirety of the parameters (the $n$ parameters of $\alpha$ and the $n^2$ of $Q$) whereas moments matching approaches only needs to match $2n - 1$ moments [BHT08] in order to obtain an estimate of the target $n+n^2$ parameters of $\alpha$ and $Q$ (based on the moments every entry of $\alpha$ and $Q$ is determined directly), hence, in this work we used moment matching methods, as they considerably reducing the computational cost.

<footnote n="3">³we assume in this work that the initial probability of the absorbing state is 0</footnote>
<figure><img src="image_47.png" /><figcaption>Figure 5.7: A degree-3 CPH: initial probability vector $\alpha$, intensities matrix $Q$, state-graph (middle) and corresponding time-to-absorption PDF (right).</figcaption></figure>

The approaches proposed in [HT07, ACH11] construct a degree-$n$ CPH distribution matching $2n - 1$ moments if the moments can be realized by a degree-$n$ CPH distribution.

Figure 5.8 shows the resulting initial probability vector $\alpha$ and transition matrix $Q$ for a degree-$n$ CPH (with $n = 1, 2$ and $3$) obtained from the Cantu hospital dataset through the BuTools package [BUT], a software tool that implements the CPH moment matching approach.

<table><caption>Figure 5.8: Results of moment matching of degree-$n$ CPH processes w.r.t. the patient's arrivals dataset.</caption><thead><tr><th>n</th><th>2n - 1<br>(moments matched)</th><th>estimated α and Q</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>α = (1)<br>Q = (-0.051453)</td></tr><tr><td>2</td><td>3</td><td>α = (0.0409229, 0.959077)<br>Q = (<table><tbody><tr><td>-0.0213213</td><td>0.0213213</td></tr><tr><td>0</td><td>-0.0570911</td></tr></tbody></table>)</td></tr><tr><td>3</td><td>5</td><td>α = (0.0284672, 0.881506, 0.0900266)<br>Q = (<table><tbody><tr><td>-0.699657</td><td>0</td><td>0</td></tr><tr><td>0.0667986</td><td>-0.0667986</td><td>0</td></tr><tr><td>0</td><td>0.0260063</td><td>-0.0260063</td></tr></tbody></table>)</td></tr></tbody></table>

From the estimated degree-$n$ CPH (Figure 5.8) we computed (through uniformisation) the PDF of the time-to-absorption for the 3 cases $n = 1, 2, 3$ and compared them with the inter-arrival PDF computed over the dataset (Figure 5.9). We point out that only the degree-3 CPH provides a good approximation of the dataset PDF, whereas the 2-phases CPH is capable of approximating only the front of the PDF while the 1-phase CPH fails to approximate the front too.

It is evident that, even if the inter-arrival time distribution is well approximated, a
<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /><figcaption>Figure 5.9: Front (left) and tail (right) of pdf of CPH distribution approximations of the inter-arrival time distribution based on moment matching.</figcaption></figure>

renewal process is significantly different from the data set under study because it cannot exhibit correlation and periodicity. Nevertheless, as we will show in Section 5.3.2, the renewal process with a proper PH distribution can give good approximation of some performance indices. In the following we propose models in which there is correlation in the inter-arrival time sequence.

### 5.2.2 Markovian arrival processes

To overcome the limitations of MRPs, i.e. the impossibility that they generate correlated arrivals, we turned to the family of continuous time Markovian arrival processes (MAP) [Neu79]. A MAP is a kind of process in which the emulation of events of interest (arrivals in our case) is governed by a background CTMC, that is, a CTMC whose states are used to model Poisson processes of given arrival intensities, and whose state-transitions allows to switch between different Poisson processes hence to model correlations. As a result a MAP allows for modelling of non-exponential and correlated inter-arrival times and may be seen as the generalization of Poisson point process.

The infinitesimal generator matrix of the underlying CTMC of an n-state MAP will be denoted by the $n \times n$ matrix $D$. Arrivals can be generated in two ways. First, in each state $i$ of the CTMC a Poisson process with intensity $\lambda_i$ is active and can give rise to arrivals during a sojourn in the state. Second, when in the CTMC a transition from state $i$ to $j$ occurs an arrival is generated with probability $p_{i,j}$. A convenient and compact notation to describe a MAP is obtained by collecting all intensities in two matrices, $D^{(0)}$ and $D^{(1)}$, in such a way that $D^{(0)}$ contains the transition intensities that do not generate an event and $D^{(1)}$ those that give rise to one, including the intensities of the Poisson processes in the diagonal of $D^{(1)}$. Moreover, the diagonal entries of $D^{(0)}$ are set in such a way that we have $D^{(0)} + D^{(1)} = D$. Accordingly, the entries of $D^{(0)}$ and $D^{(1)}$ are

$$ \forall i, j, i \neq j: \quad D_{i,j}^{(0)} = (1 - p_{i,j})D_{i,j}, \quad D_{i,j}^{(1)} = p_{i,j}D_{i,j} $$

$$ \forall i: \quad D_{i,i}^{(0)} = - \sum_{\forall j, j \neq i} D_{i,j} - \lambda_i, \quad D_{i,i}^{(1)} = \lambda_i $$
In most of the literature, the vector of steady state probabilities of the background chain, denoted by $\gamma = (\gamma_1, ..., \gamma_n)$, is used as initial probability vector of the model$^4$. In this paper we will do the same. As an example, consider the 2-state MAP described by

$$D^{(0)} = \begin{pmatrix} -3 & 1 \\ 0 & -7 \end{pmatrix} \quad D^{(1)} = \begin{pmatrix} 0 & 2 \\ 2 & 5 \end{pmatrix} \quad D = \begin{pmatrix} -3 & 3 \\ 2 & -2 \end{pmatrix} \qquad (5.1)$$

which implies that during a sojourn in state 1 no arrivals are generated while during a sojourn in state 2 a Poisson process with intensity equal to 5 is active. Moreover, a transition from state 1 to state 2 generates an arrival with probability 2/3 while transitions from state 2 to 1 are always associated with an arrival. The steady state probabilities of the background chain are $\gamma = (0.4, 0.6)$.

The well-known Poisson process is a MAP with a single state ($n = 1$). The renewal process used in Section 5.2.1 can be expressed in terms of a MAP by setting $D^{(0)} = Q$, $D^{(1)} = (-Q)\mathbb{I}\alpha$. The Markov modulated Poisson process, in which a Poisson process is active in every state of a background CTMC, is a MAP whose $D^{(1)}$ matrix contains non-zero entries only in its diagonal.

As for CPH distributions, two families of parameter estimation methods have been developed in the literature: the first based on the maximum-likelihood principle ([Ryd96, Buc03, OD09]) and the second based on matching a few statistical parameters of the arrival process. The representation given by $D^{(0)}$ and $D^{(1)}$ contains $2n^2 - n$ parameters (in every row of $D = D^{(0)} + D^{(1)}$ the sum of the entries must be zero) and it is redundant as an $n$-state MAP is determined by $n^2 + 2n - 1$ parameters ($2n - 1$ moments of the inter-arrival times and $n^2$ joint moments of consecutive inter-arrival times; see [TH07] for details). Maximum-likelihood based methods suffers from the same drawbacks described before in case of CPH distributions. In this paper we experiment with methods that belong to the second family of approaches.

A 2-state MAP is determined by 3 moments of the inter-arrival times and the lag-1 auto-correlation of the inter-arrival time sequence [BHHT08, TH07]. Our sequence has however such a lag-1 auto-correlation that cannot be realized with only 2 states. In [Hor13] a method was proposed that creates a MAP with any 3 inter-arrival time moments and lag-1 auto-correlation. This method, implemented in the BuTools package, provides the following 6-state MAP:

$$D^{(0)} = \begin{pmatrix} -0.0293 & 0.0293 & 0 & 0 & 0 & 0 \\ 0 & -0.1883 & 0 & 0 & 0 & 0 \\ 0 & 0 & -0.0103 & 0.0103 & 0 & 0 \\ 0 & 0 & 0 & -0.0989 & 0.0989 & 0 \\ 0 & 0 & 0 & 0 & -0.0989 & 0.0989 \\ 0 & 0 & 0 & 0 & 0 & -0.0989 \end{pmatrix} \qquad (5.2)$$

<footnote>$^4$Notice that $\gamma$ exists and is independent of the initial state as the background CTMC is, by definition, ergodic.</footnote>
<figure><img src="image_50.png" /><figcaption>Figure 5.10: Acf of the inter-arrival time sequence of the original data and of the MAP given in (5.2-5.3).</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 5.11: Acf of the number of arrivals per hour sequence of the original data and of the MAP given in (5.2-5.3).</figcaption></figure>

$$D^{(1)} = \begin{pmatrix} 0. & 0. & 0. & 0. & 0. & 0. \\ 0.023 & 0.1288 & 0.0002 & 0.0362 & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. \\ 0.0037 & 0.0207 & 0.0004 & 0.074 & 0. & 0. \end{pmatrix} \quad (5.3)$$

with which the steady state probability vector of the background Markov chain is $\gamma = (0.149503, 0.153343, 0.0126432, 0.22817, 0.22817, 0.22817)$.

As opposed to the approach used in Section 5.2.1, in the arrival process generated by a MAP there can be correlation between subsequent inter-arrival times. In Figure 5.10 we depict the acf in the inter-arrival time sequence generated by the above 6-state MAP and that computed on the available dataset. As guaranteed by the applied method for $n = 1$ the auto-correlation is matched exactly but then it fails to follow the auto-correlation of the data. The acf of the sequence counting the number of arrival per hour generated by the same MAP is given instead in Figure 5.11 which is very different from the that of the dataset. In the number of patients per day sequence the auto-correlation of the 6-state MAP with $n = 1$ is 0.00469433, i.e., it is negligible, as opposed to that in the data where it is 0.22.

As seen above, 6 states are necessary to match three moments and just the lag-1 autocorrelation of the inter-arrival times. This indicates that a MAP with large number of states is necessary to capture the peculiar statistical features of the patient arrival process. General purpose MAP fitting techniques are not applicable however with such large number of states.

### 5.2.3 Hidden Markov processes

A Hidden Markov model (HMM) [Rab90] can be thought of as a generalisation of a DTMC for which an external observer cannot directly see the states but only observe some output whose probability to be emitted depends on the state. In practice an HMM is characterised by: i) a set of states $S = \{s_1, ..., s_n\}$; ii) a set of possible observations $O = \{o_1, ..., o_m\}$;
iii) a $n \times n$ state-transition probability matrix $A = \{a_{ij}\}$ with $a_{ij}$ being the probability of transitioning from state $s_i$ to $s_j$; iv) a $n \times m$ state-observation probability matrix $B = \{b_{ij}\}$ where $b_{ij}$ is the probability of observing $o_j$ when the chain enters state $s_i$; v) an initial distribution over the set of states $S$ denoted by $\alpha$. Accordingly, an HMM is completely determined by the triple $(\alpha, A, B)$.

There exist three classical problems for HMM, all of which requiring a sequence of observations $O = (o_1, \dots, o_k)$. The first one is the *evaluation problem*: given an HMM by $(\alpha, A, B)$ and a sequence $O$, calculate the probability that the HMM produces $O$. The second, called *decoding problem*, consists of determining the most probable state sequence given a HMM and a sequence $O$. The third one, the *learning problem*, has only $O$ as input and is about finding such triple $(\alpha, A, B)$ with which observing $O$ is most probable.

In this paper we focus on the third type of problem and in particular we develop and study two HMMs aimed at reproducing the statistical characteristics of the patients’ arrival dataset (described in Section 5.1). We start off with a rather coarse-grained (only 3-state) but general HMM, which reveals to be capable of reproducing only marginally the auto-correlation characteristics of the patients arrivals data. Then we propose a HMM with particular underlying DTMC that shows good agreement with the dataset from a statistical point of view.

## A 3 states HMM

Here we consider a 3 state HMM in which every time slot corresponds to one hour and the possible observations are $O = \{0, 1, \dots, 14\}$ interpreted as the number of patient arrivals per hour<fnref n="5" />. This means that the sequence generated by the HMM has to be post-processed if we need to specify the exact arrival instance for each patient. This post-processing will consist of distributing the arrivals in uniform manner inside the hour.

There are 2 free parameters in $\pi$, $3 \times 2 = 6$ in $A$ and $3 \times 14 = 42$ in $B$ (because $\pi$ must be normalized and also the rows both in $A$ and $B$ must be normalized). We applied the Baum-Welch algorithm [BPSW70] to determine the optimal parameters starting from several initial parameter sets chosen randomly. With this relatively small model, the final optimal parameters obtained by the Baum-Welch algorithm are independent of the initial values (apart from permutations of the states). The obtained HMM is with $\pi = (0, 1, 0)$ and

$$ A = \begin{pmatrix} 0.865 & 0.135 & 0 \\ 0 & 0.835 & 0.165 \\ 0.134 & 0 & 0.866 \end{pmatrix} \quad B = \begin{pmatrix} 0.0007 & 0.0227 & 0.0807 & 0.1569 & 0.1756 & 0.1731 \\ 0.024 & 0.1103 & 0.2237 & 0.2249 & 0.19 & 0.1323 \\ 0.3533 & 0.3752 & 0.1989 & 0.0583 & 0.0116 & 0.0025 \\ 0.1652 & 0.0991 & 0.0625 & 0.0292 & 0.0177 & 0.009 & 0.0058 & 0.0006 & 0.0006 \\ 0.0544 & 0.0303 & 0.0052 & 0.0025 & 0 & 0 & 0.0004 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \quad (5.4) $$

<footnote n="5">⁵More than 14 patients per hour are very rare in our dataset.</footnote>
The mean number of arrivals per hour with the above HMM is 3.0862 and its variance is 5.67669 while in the data trace the mean and the variance are 3.08676 and 5.69061, respectively.

In Figure 5.12 (left) we depict the autocorrelation of the number of arrivals per hour. As one can expect, the 3-state HMM is not able to reproduce the sustained oscillation present in the data shown in Figure 5.3 (autocorrelation is negligible for $n \ge 10$) but does much better than the 6-state MAP given in (5.2-5.3) (see Figure 5.11).

## A 24 states HMM

In order to have a model that is able to exhibit oscillation in the autocorrelation function of the number of arrivals per hour, we define a 24-state HMM in which each state corresponds to an hour of the day and the transition probabilities are such that the process deterministically cycles through the 24 states. Accordingly, $A$ is 24 × 24 and its entry in position $(i, j)$ is 1 in $j = (i + 1) \pmod{24}$ and it is 0 otherwise. The initial probability vector is set to $\pi = (1, 0, ..., 0)$. As before, the possible observations are $O = \{0, 1, ..., 14\}$.

The Baum-Welch algorithm is such that parameters set to 0 initially remain 0. Consequently, the algorithm does not change the matrix $A$ and the vector $\pi$ and has an effect only on the entries of $B$. The number of parameters is larger, it is $24 \times 14 = 336$, but thanks to the deterministic behavior of the underlying DTMC the Baum-Welch algorithm determines the parameters in a single iteration. With the resulting HMM the mean number of arrivals per hour is 3.08619 while the variance is 5.67465. The hourly autocorrelation is shown in Figure 5.12 (center). This model provides very similar hourly autocorrelation to that of the data set.

We experimented also with a 24-state HMM letting the Baum-Welch algorithm to change any entry of the matrix $A$ (i.e., the initial entries of $A$ are random strictly greater than 0 and strictly smaller than 1). This way the number of free parameters is $24 \times 23 + 24 \times 14 = 888$. With this large number of parameters the Baum-Welch algorithm performs 4798 iterations and requires about 3 minutes of computation time on a standard portable computer. The resulting HMM has mean and variance equal to 3.08539 and 5.67468, respectively. The matrix $A$ has a similar structure to the one before but the probabilities of going to the next state is not 1 but a value between 0.6 and 0.99. The resulting autocorrelation structure is depicted in 5.12 (right). With the non-deterministic underlying Markov chain the autocorrelation cannot exhibit sustained oscillation and the autocorrelation vanishes after about 150 hours (i.e., 6-7 days).

## 5.3 Empirical validation of derived models

In the previous section we analysed the adequacy of each kind of *arrival* model derived from the dataset in “isolation”, that is: we used each model to generate arrivals observations and compared relevant statistics (e.g. pmf of inter-arrivals, auto-correlations of number of
<figure><img src="image_52.png" /><figcaption>Figure 5.12: Comparing the autocorrelation function (acf) of the dataset with that of the 3-state HMM (left), the 24-state HMM with deterministic (center) and non-deterministic (right) underlying Markov chain.</figcaption></figure>

arrivals per hour) computed on model's generated arrivals against that computed on the original dataset. To extend the validation of the derived models here we consider a simple ED service model (Figure 5.14) representing the services a patient may go through during his permanence at the ED and we "stimulate" it with arrivals generated on one hand by each of the arrival models derived from the dataset as well as with the arrivals of the actual dataset. We then compare meaningful service related statistics computed w.r.t. the combined arrival/service models against that computed by stimulating the service model with the original (arrivals) dataset.

For convenience we resort to the HASL model checking platform for the validation of the combined arrival/service models, hence we rely on stochastic Petri net<fnref n="6" /> encoding of both the arrival and the service process. The Petri net encoding of the Markovian arrival models derived from the dataset is pretty straightforward however the complexity of the resulting Petri net increases with that of the considered model. For the Markov renewal, the MAP and the 3-states HMM the resulting Petri net is pretty simple however the complex structure of the 24-states HMM arrival model resulted in a Petri net consisting of 63 places and 398 transitions (15 × 24 of which are immediate, 14 are timed with uniform distribution over [0, 1] and 24 are timed deterministic with delay equal to 1). Therefore, below we outline a formal Petri net encoding of periodic HMM models based on which we implemented an automatic generator of the corresponding NMGSPN models in terms of an input file for the COSMOS model checker. Similarly from the actual dataset consisting of *N* arrival events we also generate automatically its the Petri net encoding, which trivially corresponds with a Petri net consisting of *N* concurrently enabled deterministically timed transition whose deterministic delay is set to that of the corresponding arrival in the dataset.

### 5.3.1 Petri Net encoding of HMM arrival models

We introduce a Petri net formal encoding of an HMM that representing an arrival process whose arrivals are supposed to occur per discrete time unit. Such encoding allows us to

<footnote n="6">⁶specifically the NMGSPN outlined in Chapter 2.1.2.</footnote>
automatically translate an arrival HMM model into a corresponding NMGSPN that can
then be composed with a Petri net model that represents the services required by arrived
entities.

**Definition 5.2** (Petri net encoding of an *n*-states HMM). Let $\mathcal{H} = (\pi, A, B)$ be an *n*-state *m*-observations HMM whose observations per discrete time unit are assumed to follow a uniform distribution and where $\pi$ is the initial state distribution vector of $\mathcal{H}$, $A$ is the $n \times n$ state-transition probability matrix of the background chain and $B$ is the $n \times m$ state-observation probability matrix. We define $\mathcal{N}_{\mathcal{H}} = (P, T, I, O, M_0)$ the NMGSPN encoding of $\mathcal{H}$ as follows:

• $P = \{P_0, P_1, \dots, P_n, O_1, \dots, O_n, W, \text{Arrs}\}$ is the set of places whose elements are:

- $P_0$ is the initial place used to start the model according to the initial probabilities $\pi$

- $P_i$ ($1 \le i \le n$) corresponds with the *i*-th state of the background chain of $\mathcal{H}$

- $O_i$ ($1 \le i \le n$) models the triggering of observations that can occur while $\mathcal{H}$ is in the $i$-th state.

- $W$ is used to collect the observations (i.e. number of arrivals) triggered while $\mathcal{H}$ is in a given state and that has to occur in the next discrete time unit

- $\text{Arrs}$ used to gather that actual observations occurred in the last discrete time unit

• $T = T_{\mathcal{I},1} \cup T_{\mathcal{D}} \cup T_{\mathcal{I},2} \cup \{t_{\text{dist}}\}$

- $T_{\mathcal{I},1} = \{t_{0,i} : 1 \le i \le n, \pi_i > 0\}$, set of immediate transitions used to probabilistically start the process according to $\pi$, i.e. the weight of $t_{0,i}$ is given by $\pi(i)$.

- $T_{\mathcal{D}} = \{t_{i,j} : 1 \le i, j \le n, A_{i,j} > 0\}$ set of timed transition whose transitions are associated with fixed delay of one time unit and where $t_{i,j}$, if present, represents the probabilistic transition from state $i$ to $j$ as per the matrix $A$, hence $t_{i,j}$ is assigned with weight $A_{i,j}$.

- $T_{\mathcal{I},2} = \{t_{i,j}^{(o)} : 1 \le i \le n, 1 \le j \le m, B_{i,j} > 0\}$ set of immediate transitions in which $t_{i,j}^{(o)}$ models the probabilistic occurrence of $j$-th observation while in the $i$-th state of $\mathcal{H}$ as per matrix $B$, hence $t_{i,j}^{(o)}$ is assigned weight $B_{i,j}$.

- $t_{\text{dist}}$ is a timed transition associated with a Uniform distribution over time interval $[0, 1]$ used to model the actual delay of occurrence of each observation in the next time unit. Transition $t_{\text{dist}}$ is associated with infinite server policy.
* $I: T \to Bag(P)$ is the set of input arcs of transitions in T

$$I(t) = \begin{cases} \{P_0\} & \text{if } t = t_{0,i}, i \in \{1, n\} \wedge \pi(i) > 0 \\ \{P_i\}, & \text{if } t = t_{i,j}, (i, j) \in \{1, n\} \times \{1, m\} \wedge A_{i,j} > 0 \\ \{O_i\}, & \text{if } t = t_{i,j}^{(o)}, (i, j) \in \{1, n\} \times \{1, m\} \wedge B_{i,j} > 0 \end{cases}$$

* $O: T \to Bag(P)$ is the set of output arcs of transitions in T

$$O(t) = \begin{cases} \{P_i\} + \{O_i\} & \text{if } t = t_{0,i}, i \in \{1, n\} \wedge \pi(i) > 0 \\ \{P_j\} + \{O_j\}, & \text{if } t = t_{i,j}, (i, j) \in \{1, n\} \times \{1, m\} \wedge A_{i,j} > 0 \\ (j-1)\{W\}, & \text{if } t = t_{i,j}^{(o)}, (i, j) \in \{1, n\} \times \{1, m\} \wedge B_{i,j} > 0 \end{cases}$$

* $M_0 = \{P_0\} \subset Bag(P)$ is the initial marking

<figure><img src="image_53.png" /><figcaption>Figure 5.13: Petri net encoding of the 3-state HMM model given in (5.4).</figcaption></figure>

**Example 7** (Encoding of 3-states HMM derived from ED dataset). Figure 5.13 shows the NMGSPN encoding of the 3-state HMM given in (5.4) which we derived from the ED dataset. It consists of 8 places and 41 transitions, 7 of which are timed the rest immediate.
Timed deterministic transitions of $T_D$ are depicted as filled-in gray thick rectangles, the single timed uniform transition $t_{dist}$ is drawn as a black rectangle and labeled with its firing interval while immediate transitions of $T_{I,1} \cup T_{I,2}$ are depicted as thin black rectangles. The label of each timed deterministic transitions as well as that of the immediate transitions contains the “weight” of the transition (i.e. a positive real-value used to determine the probability with which the transition fires when concurrently enabled with other transitions) which corresponds with entry $A_{i,j}$, for timed deterministic transitions $t_{i,j} \in T_D$, with entry $B_{i,j}$, for the immediate “output” transitions $t_{i,j}^{(o)} \in T_{I,2}$ and with entry $\pi_i$ for initial state transitions $t_{0,i} \in T_{I,1}$, where $A, B$ are the matrices in (5.4) and $\pi = (0, 1, 0)$.

In conformance with the initial distribution of the 3-states HMM $\pi = (0, 1, 0)$, the initial marking $M_0$ has a single token in place $P_0$, and accordingly transition $t_{0,2}$ weight is set to 1. On firing of $t_{0,2}$ a token is placed in $P_2$ representing the activation of the (periodic) background chain by enabling the competing timed deterministic transitions $t_{2,3}$ and $t_{2,2}$. Notice that the timed timed transitions $t_{i,j}$ coupled with the corresponding $t_{i,i}$ ($i \in \{1, \dots, n\}, j = (i+1)\mod(n)$) model the competition between the occurrence of an observation-less transition to the “next” state of the background chain and an observation filled permanence in the $i$-th state of the chain in the next 1-time-unit, with the competition being solved probabilistically through the weights of $t_{i,j}$ and $t_{i,i}$. Whenever the competition between $t_{i,j}$ and $t_{i,i}$ is won by $t_{i,i}$ a token is added to $O_i$ activating a probabilistic choice between the immediate transitions $t_{i,k}^{(o)}$ ($k \in \{1, \dots, m\}$) which model different possible number of arrivals to occur in the next time unit, with $t_{i,0}^{(o)}$ representing 0 arrivals to be observed in the next time unit. Notice that the number of existing immediate transitions $t_{i,k}^{(o)}$ activated by a token in $O_i$ depends on the number of non-null entries in the $i$-th row of matrix $B$. For example for the 3-states HMM we know that the maximum number of arrivals per hours is $m = 14$ therefore at most 15 immediate transitions $t_{i,k}^{(o)}$ ($i.e. k \in \{0, \dots, 14\}$) can be activated by a token in $O_i$. However in the NMGSPN model in Figure 5.13 only $O_1$ is equipped with all 15 output transitions $t_{1,k}^{(o)}$ while for $O_2$ we only have 11 output transitions $t_{2,k}^{(o)}$ (e.g $t_{i,10}^{(o)}$ and $t_{i,11}^{(o)}$ are missing as the corresponding entries of $B$ are 0) and similarly for $O_3$ we only have 6 transitions $t_{3,k}^{(o)}$.

Notice that the number of occurrences of arrivals to be generated in the next time unit is modelled by setting to $k$ the multiplicity of $t_{i,k}^{(o)}$ output arc. Finally the actually occurred arrivals are represented by the tokens collected in the sink place *Arrs*.

**Petri net model of ED services.** For the sake of simplicity we assumed the ED patient flow consisting of a simple pipeline of 3 services: *triage*, *visit* and *discharge* for each of which we considered exponentially distributed service time with the following settings <fnref n="7" />: *triage*~ Exp(12), the *visit*~ Exp(6) and *discharge*~ Exp(60) (i.e. *visit* is assumed to be the slowest while *discharge* the fastest service). Figure 5.14 shows

<footnote n="7">⁷The rates have been devised from statistical analysis of the Cantù dataset.</footnote>
<figure><img src="image_54.png" /><figcaption>Figure 5.14: A simple NMGSPN model of ED services consisting of 3 pipelined services (to be composed by overlapping of place *Arrs* with an NMGSPN model of the patients arrivals).</figcaption></figure>

the NMGSPN encoding of the ED service model: notice that place *Arrs* is shared with the NMGSPN models of the patient arrivals (representing the composition of the two models). To validate each arrival model we compared their performance with that resulting by the simulation of the ED service model (Figure 5.14) with the actual arrivals of the dataset (for this we generated an NMGSPN encoding of the dataset arrivals).

### 5.3.2 KPI assessment

In order to compare the performances of the different patient arrival models coupled with the ED services we considered the following key performance indicator (KPI):

$$ \phi_1 \equiv \text{pmf of the number of patients in the ED during a periodic time-window.} $$

We formally specified $\phi_1$ by means of the Hybrid Automata Specification Language (HASL) (see Chapter 2) and assessed it through the statistical model checking platform Cosmos [BDD$^+$11a, COS].

Figure 5.15 depicts the HASL specification for $\phi_1$. The LHA consists of 4 states ($l_0, l_{on}, l_{off}, l_{end}$), 1 clock variable ($t$), 2 stopwatches ($t_a, t_p$), $M+1$ real valued variables $x_i$ ($0 \le i \le M$), whose final value is the probability that $i$ patients have been observed in the ED during the periodic time windows), plus a number of auxiliary variables and parameters.

The LHA is designed so that measuring, along a trajectory, is periodically switched ON/OFF (parameter *TP* being the ON period duration) and stops as soon as the *NP*-th period has occurred. In the initial state $l_0$ the occurrence of any transition is ignored ($l_0 \xrightarrow{ALL, t<initT, \emptyset} l_0$) for an initial transient of duration *initT* (*initT* being a parameter of the LHA) at the end of which the LHA moves to $l_{on}$ ($l_0 \xrightarrow{\sharp, t \ge initT, \emptyset} l_{on}$). In $l_{on}$ the LHA reacts to the occurrence of patients arrival as well as patients discharge events. When a patient arrives (resp. is discharged) while $i$ patients are in the ED ($l_{on} \xrightarrow{\{tarr\}, n=i \land tp \land TP \land tp \land TP \land np \land NP, n=+1, x_i+=t_a, t_a=0} l_{on}$, resp. $l_{on} \xrightarrow{\{discharge\}, n=i \land np \land NP, n=-1, x_i+=t_a, t_a=0} l_{on}$) the patients counter $n$ is incremented (resp. decremented), the duration of the last time interval on which the ED contained $i$ patients is added up to $x_i$, and the arrival stopwatch $t_a$ is reset. Any event different from a patient's arrival/discharge is ignored in $l_{on}$.
<figure><img src="image_55.png" /><figcaption>Figure 5.15: Automaton for measuring the pmf of the number of patients in the ED.</figcaption></figure>

$ (l_{on} \xrightarrow{ALL \setminus \{tarr, discharge\}} l_{on}) $. As soon as the ON-period expires ($ t_p = TP $) the LHA moves from $ l_{on} $ to $ l_{off} $ where it suspends registering the duration of different number of patients in the ED. In that respect observe that if the end of an ON-period corresponds with $ i $ patients being in the ED, the variable $ x_i $ (which accumulates the durations of $ i $ patients being in the ED) is added up with the time elapsed since the arrival of the last patient $ (l_{on} \xrightarrow{\sharp, np<NP \wedge tp=TP \wedge n=t} l_{off}) $, hence the end of the ON-period is made corresponding with the end of the duration of having $ i $ patients in the ED. In $ l_{off} $ the automata ignores any event $ (l_{off} \xrightarrow{ALL, tp<24-TP, \emptyset} l_{off}) $ while it switches back to $ l_{on} $ as soon as the reciprocal (w.r.t. 24 hours) of the ON-period duration elapsed $ (l_{off} \xrightarrow{\sharp, tp=24-TP, \{t_a=0\}} l_{on}) $ and in so doing the timer $ t_a $ (which stores the duration of the occupation at $ i $ patients since the last arrival/departure) is reset so that it can correctly be used in the freshly started ON-period. Finally, the LHA stops monitoring as soon as NP ON-periods have been observed along the monitored trajectory: at that moment each $ x_i $ is normalised w.r.t. the sum of all $ x_i $, hence on ending the monitoring of trajectory $ x_i $ is assigned with the probability that the $ i $ patients have been observed in the ED during NP observed periods $ (l_{on} \xrightarrow{\sharp, np=NP, \{x_0=\frac{x_0}{\sum_{i=0}^{M} x_i}, \dots, x_M=\frac{x_M}{\sum_{i=0}^{M} x_i}\}} l_{end}) $. Observe that such an LHA can also be used for “non-periodic” measures: it suffices to set $ TP = 0 $. The HASL specification for $ \phi_1 $ is completed by the list of HASL expressions $ AVG(last(x_i)) $ which indicate that a confidence interval for average of the last value that $ x_i $ has at the end of an accepted trajectory is computed by COSMOS.

Figure 5.16 depicts plots computed with the COSMOS tool and resulting from assessing specification $ \phi_1 $ (i.e., the pmf of the number of patients in the ED) against the NMGSPN models of the patient arrivals coupled with the 3-services model of the ED (Figure 5.14).
<figure><img src="image_56.png" /><figcaption>Figure 5.16: Pmf of the num. of patients in the ED computed through 3-phases PH renewal model, 6-state MAP, 3-state HMM and 24-state HMM versus dataset model over the whole day (top left), during a low-arrival hour (top right) and high-arrivals hour (bottom).</figcaption></figure>

Plotted results refer to a 1-year (365 days) observation window and have been computed as 99.99% confidence interval of $10^{-3}$ width. The top-left plot refers to the pmf measured without taking into account any specific time-window over the day (continuous measuring over 24h for 365 days), while the top-right and bottom plots refer to periodic measuring over a 1-hour period, at low arrival intensity (from 2am to 3am, top-right plot) and at high arrival intensity (from 11am to 12pm, bottom plot), respectively. The obtained results witness the clear advantage of the HMM24 model over the rest: if when no specific hour of the day is considered the pmf of the 3-phases renewal model and of the 6-states MAP provide an acceptable approximation of the pmf computed w.r.t. to dataset (red dashed plot) and the pmf of the HMM3 and HMM24 are essentially indistinguishable, this is no longer the case when the pmf is computed on a specific time window as shown in the top-right and bottom plots. Only the pmf of the HMM24 matches that of the dataset during a low-arrivals [2am,3am] window and a high-arrival [11am,12pm] window: the pmf of the renewal and map models instead are essentially identical to those measured continuously while that for the HMM3 exhibits a slight tendency towards matching the pmf of the dataset at low-arrivals but completely fails to do so when arrivals become intense.
## 5.4 Perspectives

The contribution presented in this chapter tackles the problem of deriving stochastic models capable of reproducing the periodic character of observations of a real-life phenomenon exhibiting periodicity. Based on a specific dataset that registered observations of patients arrivals to the emergency department of an hospital in northern Italy we derived different instances of Markovian models which we proved able, to different extent, to generate arrival events reproducing relevant statistical indicators assessed on the dataset itself. More specifically we showed that, amongst different families of Markov models, the class of doubly stochastic processes such as, e.g. Hidden Markov models (HMMs) is worth considering as it has the potential to accurately reproduce the periodic character of the dataset, given it is “properly” set. In fact, generally speaking, classical hidden Markov models fail to adequately capture periodic behaviours because, even if a large number of states is used, after an initial transient phase where the process might be able to mimic the periodicity of the observed phenomenon, as the behavior of the Markov chain, at steady state, “smooths out” and the correlation of the probability of a given state and the elapsed time disappears. To face this problem we forced the background Markov chain of an HMM to be periodic, i.e., any return to any state occurs in multiples of a given number of time steps. Such chains are capable to act as a clock, and hence maintain the periodic pattern of the phenomenon we aim to study, as we demonstrated for the patients arrival dataset. Although promising, the approach outlined in this chapter is somehow limited as it only considers arrival of discrete quantities (arrivals of patients to an ED, or cars to a traffic light), although many real life systems are rather concerned with continuous quantities observations (e.g. energy produced by a photovoltaic system or the amount of yogurt produced in a dairy industry). Therefore it’d be worth to consider how to generalise the proposed approach so that it can cope with other family of periodic phenomena. In practice that implies studying what are the implications of considering different kind of arrival processes while sojourning in the states of the periodic background chain of an HMM. A number of different situations seem worth considering in this respect, including using of deterministic arrival intensity as opposed to random arrival intensity and also using of continuous and/or discrete phase-type distributions to models observations of continuous or discrete observations.

In particular a number of open questions need to be tackled related to parameter estimation algorithms for HMMs with periodic background chain, including: 1) would it be possible to adapt/improve expectation-maximization parameter estimation approaches to HMMs whose background chain is periodic? (experience suggests that special background structure, like a periodic chain, can lead to more efficient calculation for what concerns the expectation step); 2) would it be possible to adapt/improve moment fitting approaches to HMM when phase-type distributions are used to generate observations while sojourning in the states of the periodic background chain? Those research questions have been integrated within a PhD project proposal that Andras Horváth and myself have written and for which we are currently looking for suitable PhD candidate.
# Chapter 6

## Other contributions

In this chapter we provide the reader with a brief overview of other contributions not included in this manuscript.

**Compositional approaches for stochastic models.** Together with Andras Horváth from University of Turin (Italie) we worked on compositional approaches for assessing probabilistic properties of structured Markov chain models. This entails looking at how to exploit the compositional structure of a model in order to alleviate the memory needs of techniques for assessing the model's properties. In a preliminary work [BH08] we focused on the *Boucherie* product processes, a specific class of *product form* continuous time Markov chains and showed that the compositional constraints that lead to the product form result for such class, can be exploited, to a limited extent, in the model-checking problem as well, leading to a decomposed semantics for a fragment of the Continuous Stochastic Logic. As a follow up we looked at a more general class of structured Markov chains those whose state-space can be partitioned in an arbitrary number of disjoint subsets, called macro-states, where macro-states represent the concurrent execution of independent jobs. These kind of chains typical arise from commonly used formalisms such as stochastic Petri nets or stochastic process algebras where the stochastic activities are modelled through phase type distributions. We showed that the computation of the probability that this kind of process passes through a given sequence of macro-states, i.e. the probability of a path consisting of macro-states, can be calculated in a memory efficient manner by means of Laplace transform techniques [BH09a, BH09b].

**Transient analysis of stochastic timed automata networks.** Together with Nathalie Bertrand of INRIA Bretagne/Atlantique (France), Enrico Vicario and Marco Paolieri from University of Florence (Italie) and Andras Horváth from University of Turin (Italie) we introduced a novel formalism for modelling of concurrent stochastic systems, namely that of interacting stochastic timed automata (STA) and equipped it with an operational
semantics that allow for transient-state analysis [BBH$^{+}$13]. We introduced STA as an extension of timed automata in which both delays and actions are chosen at random through probability distributions. By imposing specific rules on the transition firing, i.e. from the current state first we randomly select the next action to occur and then we sample the sojourn time (i.e. the delay of occurrence of the selected transition) we ensure that the underlying semantics can be described by means of the stochastic state classes method [BPSV05], that is, a method formerly introduced in the context of stochastic timed Petri nets and that allows for evaluation of steady-state, transient-state and quantitative model checking queries of stochastic processes with regeneration points.

**Improved Bayesian parameter estimation.** Chapter 3 reports about merging of Bayesian parameter estimation with statistical model checking yielding a novel Approximate Bayesian Computation (ABC) algorithm for parametric verification of temporal logic specifications. That work actually originated from a previous project in which in collaboration with Konstantinos Koutroumpas, Paul-Henry Cournède and Irene Votsi we developed a novel variant of the ABC sequential Monte Carlo (ABC-SMC) scheme, a popular algorithm for data-based parameter estimation of a mathematical (deterministic or stochastic) model. Given a (set of) observation(s) $y_{obs}$ of a given parametric system $\mathcal{M}_{\theta}$ the goal of ABC-SMC approaches is to estimate the parameters $\theta$ that yields a good approximation of the observations. ABC-SMC schemes operate iteratively by progressively increasing the precision of the estimate and employ Kernel distributions to propagate (i.e. “move’ around”) the parameters within the parameter space. The propagation phase of ABC-SMC turns out to be crucial to the performance of the algorithm, and in this respect the chosen sequence of Kernel distributions can be fundamental. In our contribution, which has been published on the Bioinformatics journal [KBVC16], we introduced a novel version of the ABC-SMC scheme which employs a specific family of stochastic processes, namely the Dirichlet process mixture, as Kernel distributions to move the parameters during the parameters space exploration. We showed that this results in a reduced runtime (i.e. a faster convergence) of ABC-SMC scheme particularly when the posterior distribution of the parameters to estimate is multi-modal. We validated the proposed approach w.r.t. a complex biological network model, that of the Wnt signalling pathway.

**Temporal logic based reverse-engineering of genetic regulatory networks.** As part of the PhD thesis work of Emmanuelle Gallet, which I co-supervised with Pascale Le Gall, we developed a method for reverse-engineering of parametric, discrete-state, genetic regulatory network (GRN) models. Generally speaking GRN models are sensed to capture the possible evolutions of a population of interdependent proteins whose expression is regulated by *activation/inhibition* of the corresponding genes. In their discrete abstraction a GRN is modelled by a *transition system* (also called a *dynamics*) whose states encodes the proteins’ *counts* (or *levels*) and whose transitions represent the effect that
concurrent regulation dependencies have on the proteins' counts. In this respect the René Thomas formalism constitutes an underspecified (parametric) form of GRN model, one in which level dependent genes' inter-dependencies are represented through an interaction graph without specifying what is the combined effect that (potentially concurrent) regulators have on the overall state of the network. Such an effect depends on the so-called parameters of the Thomas' GRN, which, once instantiated, allow to obtain a corresponding *dynamics* (i.e. transition system). A relevant problem considered in the literature is that of the identification of GRNs' parameters instances yielding a transition system that fulfil specific behaviours expressed by formulae of the LTL logic. Such a problem suffers from combinatorial explosion as the number of parameters is exponential in the number of genes (and levels) of the network. Solutions have been proposed (e.g. Bernot et al. [BCRG04]) that exhaustively search the parameter space and that are implemented by the SMBionet tool [KCAB09], however they are affected by scalability issues. Within Emmanuelle Gallet's PhD thesis we proposed an alternative approach in which we employ *symbolic execution* techniques to solve a parametric LTL model checking problem, the one resulting by considering the product of a parametric (Thomas) GRN with the Büchi automata corresponding to an LTL formula. By so doing we avoid the costly parameter space exhaustive search and obtain a tool that scale up remarkably better than SMBionet. Results of this research line have lead to the following publications [GMGB14, GMGB15], other than to realisation of a prototype software tool named SPuTNIk.

**Framework for formal modelling of manufacturing systems.** The design of modern industrial production systems is strongly affected by product quality and delivery reliability requirements: the ability to guarantee that products are issued within given time deadlines and that they match given quality standards are essential factors. As part of a research project in a recent collaboration with Andras Horváth (University of Turin) we developed a probabilistic model checking framework for formal modelling and performance analysis of so-called *synchronous production lines* that is, production systems consisting of a multi-stage working line where a workpiece is progressively transformed into the end product through a sequence of processing stages performed by dedicated pipelined workstations (i.e. a workpiece outputted by one workstation is fed in to the next downstream station). Workstations are assumed to be unreliable (can break down hence requiring repairing) and are separated by finite-size buffers (where pieces are temporally stored throughout the intermediate phases of manufacturing) which may induce the blocking of a workstation (e.g if the downstream buffer is full). Furthermore the overall dynamics of such systems is assumed to be strongly synchronous: the state of machines (working or broken) as well as that of the buffers evolve synchronously according to the elapsing of discrete units of time.

It has been shown [CT12, CAH14, AHC14, CHA15, ACHG16] that this kind of production systems are suitable to be modelled through discrete-time Markov chains (DTMCs) and that relevant key performance indicators (KPIs), such as the *lead time distribution* (i.e. the probability distribution of the time taken by the system to output one product),
can be assessed by application of classical transient-analysis of Markov chain models. The downside of those kind of approaches is that they are based on an explicit matrix representation of the DTMC that represents the considered production line, hence they provide little support from a modeller point of view. Furthermore the type of KPIs that one can take into account with this approaches is somehow limited.

In order to tackle these limitations we introduced (two versions of) a model-checking based framework which allows one to automatically generate a DTMC model, through an adequate modelling formalism, for a production line of arbitrary size while, at the same time, a number of relevant KPIs expressed in terms of some temporal logic formalism are also generated. The generated model and KPIs can then be promptly used to run model checking experiments relying on the model checking platform targeted by the framework hence proving an effective means for customised performance analysis of production lines. A first version of such framework [BH21a] is aimed at the PRISM model checking platform and relies on the Continuous Stochastic Logic for expressing relevant KPIs. With this we provide the user with a generator of PRISM models as well as of relevant KPIs expressed as CSL formulae and we further take advantage of model checking to validate the generated models through the verification of a number of sanity checks properties. Although effective this first version of the production line modelling framework is undermined by the limited expressiveness of the CSL property language which restrain the KPIs that one can take into account. To overcome this issue we developed a second version of such framework [BH21b] which targets the COSMOS model checking platform hence taking advantage of the much more expressive HASL [BBD$^+$15b] property language. Since COSMOS models are expressed through a stochastic extension of the Petri nets formalism the model generator of this version of the framework generates stochastic Petri nets encoding of production lines of arbitrary size as well as a number of sophisticated KPIs encoded as HASL formulae.

**Formal analysis of wireless network protocols.** The exchange of information between computing nodes over a wireless network is regulated through Medium Access Control (MAC) algorithms implementing so-called carrier sense multiple access collision avoidance (CSMA-CA) schemes. The basic idea is to employ probabilistic *backoff* procedures in order to reduce collisions between concurrently transmitting nodes, therefore letting a transmitting node *hold on*, for a probabilistically chosen duration, if the shared (wireless) medium is already occupied. This backoff mechanism is at the basis of popular IEEE wireless network standards, ranging from the basic 802.11 standard that addresses wireless local area networks (WLAN), to its *prioritised* extensions in which prioritised traffic classes are introduced to support QoS in broadband wireless network, as is the case with the 802.16 (WiMAX) standard or with the 802.11p standard for Vehicular Ad Hoc Networks (VANETs), i.e. a relevant class of networks in the realm of Intelligent Transportation Systems.

Performance analysis of wireless network systems based on IEEE MAC protocols at-
tracted the attention of researchers also in the formal methods community. Inspired from a study by Jeremy Sproston et al. [KNS02] where performance properties of the basic access (BA) version of the 802.11 MAC were formally assessed through the PRISM model checker over a (toy) single-hop wireless network, together with Alice Miller (University of Glasgow) we addressed formal performance analysis of the more complex ready-to-send/clear-to-receive (RTS/CTS) version of 802.11 MAC and considered the comparison against its energy saving version, namely S-MAC i.e. a backoff collision avoidance mechanism operating on top of a sleep-awake cycle through which the nodes’ radio is periodically turned off so to preserve the battery consumption in energy constrained wireless sensor nodes. To account for more realistic situations we considered a larger three-hops network and assessed how performance indicators (transmission latency) trade off against energy indicators (energy consumption) for transmission of data packets over a 3-hops [BM06]. The main limitation of such performance modeling effort is in terms of scaling as the complexity of the underling discrete-time Markov chain model (encoded in the PRISM model checking language) is affected by the network size yielding the necessity to generate very large, error prone, PRISM code which also quickly lead to saturating the memory when model checking algorithms are run on them.

In order to allow the modellers to take into account for larger networks consisting of an arbitrary number of stations, in collaboration with Nicolas Vasselin and Benoit Barbot (Université Paris-Est Créteil), we developed a novel modelling framework for performance analysis of protocols of the IEEE 801.11 MAC family and specifically suitable for comparing the standard RTS/CTS 802.11 MAC for wireless local area networks (WLAN) against its prioritised version used in Vehicular Ad Hoc Networks (VANETs) [BBV19]. This framework relies on statistical model checking (the COSMOS platform) hence removes the memory limitations numerical model checking approaches suffer of, and is based on a high-level stochastic Petri net formalisms, namely the so-called symmetric stochastic nets, which provides for a compact representation of complex systems consisting of a multitude of components copies. Based on the expressiveness of coloured Petri nets formalism the framework's is highly configurable as, for example, the number of stations the network consists of corresponds with a *color domain* of the Petri net model. Therefore the framework can account for i) an arbitrarily large number of stations, ii) different traffic conditions (saturated/non-saturated), iii) different hypothesis concerning the shared channel (ideal/non-ideal).

Another relevant class of protocols is the 802.16 WiMAX standard, which addresses QoS needs in broadband networks by means of dedicated mechanisms, such as the Ad-mission Control (AC), the Traffic Policer (TP) and the Traffic Shaper (TS) with different responsibility in the QoS management process. The AC which resides on the base-station (BS), establishes whether to accept or reject an incoming connection request according to the negotiated QoS parameters (maximum sustained traffic rate, minimum reserved traffic rate, and maximum latency), the TP, which resides on subscriber-stations (SSs), drops incoming requests that violate a pre-negotiated QoS agreements, while the TS, as opposed
to the TP, delays (rather than dropping) incoming requests violating a pre-negotiated QoS agreement. More specifically the role of TS is to send packets into the network according to the QoS needed by each service class and avoiding congestion. In collaboration with Lynda Mokdad (University Paris-Est Créteil) and Jalel Ben-Othman (University Paris 13) we analysed QoS tradeoffs between strict (deterministic) prioritised TS versus a novel probabilistic TS policy [BBM12].

Another relevant issue related to the security of (the energy constrained) wireless sensor networks is that of detection of denial-of-service (DoS) attacks. A DoS attack takes place when an intruder manages to get control of a network node and inject a mischievous behaviour for example by flooding the network with unwanted traffic leading to draining of the nodes’ batteries hence ceasing the network functioning. Mechanisms for detection of DoS attacks in hierarchically clustered WSNs assign a node a *control* responsibility of analysing the traffic inside a cluster and to send warnings to the cluster head whenever an abnormal behavior (e.g., high packets throughput) is detected. It is therefore paramount to avoid control nodes (cNodes) to run out of battery leaving the network unprotected. As part of the PhD thesis work of Quentin Monet (University Paris-Est Créteil), we proposed new mechanisms for DoS detection. In a first work [BMM13] we proposed a novel dynamic cNodes displacement scheme according to which cNodes are periodically elected, at random, among ordinary nodes of each cluster in order to improve the cNode battery lifetime hence the network’s w.r.t. the classical mechanisms whereby cNodes are statically displaced in strategic positions within the network topology. In a follow up work we extended that approach by proposing additional dynamic cNode election schemes in which residual energy level is also taken into account amongst the criteria for periodically electing a cluster’s cNode. We demonstrated the effectiveness of the newly introduced dynamic cNode schemes through experiments run on the network simulation platform ns-2 whose outcomes showed improvements of the load balancing in the network, while maintaining good detection coverage, w.r.t. static cNode schemes [MMB$^{+}$17].

**Jump-diffusion approximated analysis of density dependent Markov chains.**

*Density dependent (continuous time) Markov chains* (DDMC) is a subclass of CTMC used to model populations of interacting objects whereby the intensity of the interactions (state-transitions) can be expressed as a function of a parameter *N* that represents the density of the objects present in a given volume (e.g a cell’s volume), which makes chemical reaction networks (i.e. systems consisting of biochemical species interacting within a cell of given volume) a typical case of system suitable for DDMC modelling. It has been shown that DDMCs can be approximated by less computationally hungry processes with a “fluid” state space, i.e. either by ordinary differential equations (ODEs) or by stochastic differential equations (SDE) [Kur76]. Differently from ODEs approximation, SDE approximation of DDMC, known as *diffusion approximation*, preserve the stochasticity of the original DDMC, which on the contrary is lost with ODE approximation. A limitation of
the DDMC diffusion approximation is that it can be applied only to models where the probability of reaching a boundary of the state space of the process is negligible, which motivated the introduction of so-called *jump diffusion* approximation where essentially the SDE used to approximate a given DDMC is enriched with a (discrete) Poisson counting component that gives rise to jumps that mimic the behavior of the original CTMC at the boundaries [BS17, BBH$^{+14}$].

In collaboration with Marco Beccuti, Enrico Bibbona, András Horváth, Roberta Sirovich and Jeremy Sproston from University of Turin we developed a framework to estimate the probability of timed properties via *jump diffusion approximation* of large DDMC models of biological pathways [BBB$^{+17}$]. More specifically we demonstrated that statistical estimation of probabilistic time-bounded properties against DDMC models of chemical reaction networks can take advantage of jump-diffusion approximations. Specifically we showed that by replacing (costly) “exact” trajectories sampled from the DDMC through Gillespie’s stochastic simulation algorithm (SSA) with “approximated” trajectories sampled through jump-diffusion approximation one can obtain a remarkable speed up in estimating the probability of time bounded measures.
# Chapter 7

## Conclusion and perspectives

This thesis manuscript summarised some of my contributions in the domain of formal modelling for complex systems featuring probabilities and parameters. More specifically Chapter 2 gave an overview of an expressive, hybrid automata based, framework for statistical model checking that targets a very generic class of stochastic processes, while Chapter 3 tackled the problem of reverse-engineering of parametric Markov chains model, by introducing a framework that allows the modeller to tune a stochastic model's parameters so that the model's meets given, formally expressed, requirements. Chapter 4 then showed the relevance of the approaches discussed in Chapter 2 and Chapter 3, by considering a non-trivial problem, that is, the formal analysis of stochastic oscillators. Finally Chapter 5 presented a recent contribution tackling the problem of deriving stochastic models capable of reproducing the periodic behaviour of a dataset representing an observed phenomenon. Ongoing work and perspectives were given at the end of each chapter, I recall the main ones here and more generally discuss directions I envision for my research in the coming years.

**Speeding up statistical model checking.** Despite low *space complexity* SMC runtime can be remarkably high due to the cost of sampling (large amount of) paths through stochastic simulation algorithms. Indeed existing SMC frameworks use *exact* stochastic simulation sampling schemes, such as Gillespie's [Gil77], through which each *jump* in a path as a random duration and is obtained through a *stochastic race* between all enabled events. In this manner the unfolding of a path can become extremely costly particularly for models with a high degree of concurrency. In an attempt to reduce the simulation runtime *approximate* stochastic simulation algorithms, the basic idea being to speed up the simulation process by replacing the “exact” random long jumps of a path with *discrete jumps* of arbitrary fixed duration yielding an approximated path. Such approximated algorithms have been introduced in the system's biology community as means to simulate large chemical reaction networks but, at the best of my knowledge, have not been considered in the SMC context. I have been starting looking at the problem of taking advantage of
approximated simulation scheme within SMC, the main problem to tackle in this respect is to figure out how to compensate or limit the inherent error that approximated paths induce on the estimation of the satisfaction probability of a temporal logic formula. A promising approach in this respect is the multi-level simulation algorithm [LBGY16] which based on the *thickening property* of Poisson random variables, allows for cleverly generating, in parallel, paths at different level of precision and progressively compensating the error of approximated paths hence to obtain precise estimations at lower computational cost than exact schemes.

**Evolution of the automaton-ABC parametric verification framework.** The parametric model checking framework presented in Chapter 3 is, in its current formulation, limited to non-nested spatio-temporal properties. Even though this kind of properties is that found in most case studies in the literature it would be worth to work out a generalisation to the full MITL spectrum, similarly to what has been done in studies that address the robustness of MITL against real valued signals of non-probabilistic models [DM10]. This would entail to investigate the possibility to come up with a satisfiability distance definition based on the recursive syntax of MITL formulae, and then to work out a machinery for assessing distances of nested formulae against trajectories sampled from a CTMC model, which might not be an easy task considering that in automata-based model checking approaches, the characterisation of automata capable of accepting paths that satisfy nested temporal formulae turned out not an easy task. Another relevant evolution I wish to engage with is that of integration of the automaton-ABC framework within the COSMOS tool. To this day only a prototype implementation (written in the Julia programming language https://gitlab-research.centralesupelec.fr/2017bentrioum/tcs2021) of the parametric verification framework exists, which although computationally efficient is little user friendly (i.e. the necessary *distance automata* must be user as objects of some Julia class). The integration within the COSMOS platform is highly desirable as it would allow the user to rely on the HASL formalism all together, hence allowing one to express models as stochastic Petri nets and to encode the distance automata directly as a LHA.

**Stochastic process mining.** As part of the PhD work of Pierre Cry (first year PhD student at MICS laboratory of CentraleSupélec I am co-supervising since October 2022) I recently got interested in the extensions of process mining to the realm of stochastic models. Process mining algorithms are concerned with the derivation of formal models, usually in form of a workflow Petri net, capable of matching observations stored in *event logs* (collections of traces consisting of sequences of events witnessing executions of the observed process). The quality of mined models is established by comparing the behavior seen in the event-log with that exhibited by the model, w.r.t. different criteria including *fitness* (how much of the event-log behavior is reproduced by the model) and *precision* (how much of the model's behavior is not contained in the event-log). The vast majority of process mining frameworks aim untimed, non-probabilistic behaviour, i.e., they only account for the temporal order of events observed in the log while completely ignoring
timing and probability of event occurrence. The extension of process mining to stochastic models seems little explored so far in the literature. From a theoretical point of view two type of problems seem relevant for extending process mining to the stochastic realm. First one have to adapt mining algorithms so that they issue stochastic models based on “timed” event logs, something which, intuitively, seem to require a fitting phase where the stochastic distributions to be used in the stochastic (Petri net) model are established based on the timing information contained in the log. Secondly conformance checking of mined models will also need to be adapted to the stochastic case, and that likely will entail adaptation of methods for comparing stochastic languages, such as, e.g. the (Wasserstein distance) earth movers method.

**Data-driven predictive modeling of periodic phenomena.** In Chapter 5 we presented preliminary results about the derivation of models that reproduce the periodicity exhibited by observations stored in a dataset. Following a rather empirical approach we showed that the periodicity of the data can be accurately mimicked by imposing a corresponding “periodicity” in the structure of the background chain of a hidden Markov model and more specifically that sustained periodicity can be obtained by means of a “clock” background chain, i.e., a “deterministic” chain consisting of as many states as the number of (discretised) intervals in which observations are grouped and whose states form a loop through a sequence of probability 1 transitions. These preliminary findings calls for further developments. First of all it is worth investigating whether the use of less constrained periodic background chain (not necessarily clock chains) can still lead to reproduce the sustained periodicity of the data. Also through our experiments we observed that parameter estimation (expectation maximisation) algorithms used to fit the HMM to the dataset seem to converge immediately, which suggests that an explicit solution for parameter estimation might exist for this specific class of models, that is, for HMM with periodic background chain. Furthermore the overall approach should be generalised so to address other cases, such as, that of continuous data as opposed to discrete data (i.e. number of arrivals) we took into account so far. In this respect it will also worth considering the use of more generic stochastic processes in order to generate the actual observations in each state of the background chain. So far, since we addressed discrete dataset only, we only considered point Poisson process as the type of process to generate the observations. It’d be worth considering Phase type distribution both discrete (in case of discrete dataset) and continuous (in case of continuous dataset). That would entail also investigating how known expectation maximisation and moments fitting methods would cope with fitting the parameters in the case of this particular class of (periodic) structured HMM. In collaboration with András Horváth (University of Turin) we have submitted a PhD project proposal to a funding program for a co-supervising (cotutelle) PhD thesis between University of Turin and University Paris Saclay for which we hope to get funded.

The expected outcome for such a PhD project is twofold. From a practical point of view we expect that an effective framework is developed for mining of models based on periodic
dataset. Such framework should allow for (semi) automatic derivation of stochastic models
that accurately reproduce the periodic character of the dataset and should be validated
through different datasets we have access to, both from hospital emergency department
as well as from manufacturing industry domain. From a theoretical standpoint we expect
that the working out of such a framework will likely lead to novel theoretical results, such
as, for example, the introduction of novel parameter estimation procedure optimised to the
case of structured periodic stochastic models.
Bibliography

[ABB<sup>+</sup>13] Elvio Gilberto Amparore, Paolo Ballarini, Marco Beccuti, Susanna Donatelli, and Giuliana Franceschinis. Expressing and computing passage time measures of GSPN models with HASL. In José Manuel Colom and Jörg Desel, editors, *Application and Theory of Petri Nets and Concurrency - 34th International Conference, PETRI NETS 2013, Milan, Italy, June 24-28, 2013. Proceedings*, volume 7927 of *Lecture Notes in Computer Science*, pages 110–129. Springer, 2013.

[ABDF11] E.G. Amparore, M. Beccuti, S. Donatelli, and G. Franceschinis. Probe automata for passage time specification. In *Quantitative Evaluation of Systems (QEST), 2011 Eighth International Conference on*, pages 101–110, sept. 2011.

[ACD91] R. Alur, C. Courcoubetis, and D. Dill. Model-checking for probabilistic real-time systems. In *ICALP'91*, LNCS 510, 1991.

[ACH11] F. Avram, D.F. Chedom, and A. Horváth. On moments based Padé approximations of ruin probabilities. *Journal of Computational and Applied Mathematics*, 235(10):3215–3228, 2011.

[ACHG16] Alessio Angius, Marcello Colledani, András Horváth, and Stanley B Gershwin. Analysis of the lead time distribution in closed loop manufacturing systems. *IFAC-PapersOnLine*, 49(12):307–312, 2016.

[ACHH92] R. Alur, C. Courcoubetis, T. A. Henzinger, and P.-H. Ho. Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems. In *Hybrid Systems*, LNCS 736, 1992.

[ACHH93] Rajeev Alur, Costas Courcoubetis, Thomas A. Henzinger, and Pei Hsin Ho. Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems. In Robert L. Grossman, Anil Nerode, Anders P. Ravn, and Hans Rischel, editors, *Hybrid Systems*, pages 209–229, Berlin, Heidelberg, 1993. Springer Berlin Heidelberg.
[AD94] R. Alur and D. L. Dill. A theory of timed automata. *Theoretical Computer Science*, 126(2), 1994.

[AHC14] A. Angius, A. Horváth, and M. Colledani. Moments of accumulated reward and completion time in Markovian models with application to unreliable manufacturing systems. *Performance Evaluation*, 75:69–88, 2014.

[AMBC$^{+95}$] M. Ajmone Marsan, G. Balbo, G. Conte, S. Donatelli, and G. Franceschinis. *Modelling with Generalized Stochastic Petri Nets*. John Wiley & Sons, 1995.

[ANJB10] Matthew A Nunes and David J Balding. On optimal selection of summary statistics for approximate bayesian computation. *Statistical applications in genetics and molecular biology*, 9:Article34, 01 2010.

[AP18] Gul Agha and Karl Palmskog. A survey of statistical model checking. *ACM Transactions on Modeling and Computer Simulation*, 28:1–39, 01 2018.

[ASSB96] Adnan Aziz, Kumud Sanwal, Vigyan Singhal, and Robert Brayton. Verifying continuous time markov chains. In Rajeev Alur and Thomas A. Henzinger, editors, *Computer Aided Verification*, pages 269–276, Berlin, Heidelberg, 1996. Springer Berlin Heidelberg.

[ASSB00] A. Aziz, K. Sanwal, V. Singhal, and R. Brayton. Model-checking CTMCs. *ACM Trans. on Computational Logic*, 1(1), 2000.

[Bal15a] Paolo Ballarini. Analysing oscillatory trends of discrete-state stochastic processes through hasl statistical model checking. *International Journal on Software Tools for Technology Transfer*, pages 1–22, 2015.

[Bal15b] Paolo Ballarini. Analysing oscillatory trends of discrete-state stochastic processes through HASL statistical model checking. *STTT*, 17(4):505–526, 2015.

[BB22] Paolo Ballarini and Benoît Barbot. Cosmos: Evolution of a statistical model checking platform. *SIGMETRICS Perform. Evaluation Rev.*, 49(4):65–69, 2022.

[BBB$^{+17}$] Paolo Ballarini, Marco Beccuti, Enrico Bibbona, András Horváth, Roberta Sirovich, and Jeremy Sproston. Analysis of timed properties using the jump-diffusion approximation. In Philipp Reinecke and Antinisca Di Marco, editors, *Computer Performance Engineering - 14th European Workshop, EPEW 2017, Berlin, Germany, September 7-8, 2017, Proceedings*, volume 10497 of Lecture Notes in Computer Science, pages 69–84. Springer, 2017.
[BBC19] Mahmoud Bentriou, Paolo Ballarini, and Paul-Henry Cournède. Reachability design through approximate bayesian computation. In *International Conference on Computational Methods in Systems Biology*, pages 207-223. Springer, 2019.

[BBC21] Mahmoud Bentriou, Paolo Ballarini, and Paul-Henry Cournède. Automaton-abc: A statistical method to estimate the probability of spatio-temporal properties for parametric markov population models. *Theor. Comput. Sci.*, 893:191-219, 2021.

[BBD$^{+15}$a] P. Ballarini, B. Barbot, M. Duflot, S. Haddad, and N. Pekergin. Hasl: A new approach for performance evaluation and model checking from concepts to experimentation. *Performance Evaluation*, 90:53 – 77, 2015.

[BBD$^{+15}$b] Paolo Ballarini, Benoît Barbot, Marie Duflot, Serge Haddad, and Nihal Pekergin. HASL: A new approach for performance evaluation and model checking from concepts to experimentation. *Perform. Eval.*, 90:53-77, 2015.

[BBD$^{+15}$c] Paolo Ballarini, Benoît Barbot, Marie Duflot, Serge Haddad, and Nihal Pekergin. Hasl: A new approach for performance evaluation and model checking from concepts to experimentation. *Perform. Eval.*, 90:53-77, 2015.

[BBDH17] Benoît Barbot, Béatrice Bérard, Yann Duplouy, and Serge Haddad. Statistical Model-Checking for Autonomous Vehicle Safety Validation. In *Conference SIA Simulation Numérique*, Montigny-le-Bretonneux, France, March 2017. Société des Ingénieurs de l'Automobile.

[BBH$^{+13}$] Paolo Ballarini, Nathalie Bertrand, András Horváth, Marco Paolieri, and Enrico Vicario. Transient analysis of networks of stochastic timed automata using stochastic state classes. In Kaustubh R. Joshi, Markus Siegle, Mariëlle Stoelinga, and Pedro R. D'Argenio, editors, *Quantitative Evaluation of Systems - 10th International Conference, QEST 2013, Buenos Aires, Argentina, August 27-30, 2013. Proceedings*, volume 8054 of *Lecture Notes in Computer Science*, pages 355-371. Springer, 2013.

[BBH$^{+14}$] Marco Beccuti, Enrico Bibbona, András Horváth, Roberta Sirovich, Alessio Angius, and Gianfranco Balbo. Analysis of Petri net models through stochastic differential equations. In *Proc. PETRI NETS 2014*, volume 8489 of LNCS, pages 273-293. Springer, 2014.

[BBM12] Paolo Ballarini, Jalel Ben-Othman, and Lynda Mokdad. Quantitative verification of wimax traffic shaping solutions. In Elhadi M. Shakshuki and Muhammad Younas, editors, *Proceedings of the 3rd International Conference*
on *Ambient Systems, Networks and Technologies (ANT 2012)*, the 9th *International Conference on Mobile Web Information Systems (MobiWIS-2012)*, *Niagara Falls, Ontario, Canada, August 27-29, 2012*, volume 10 of *Procedia Computer Science*, pages 1026-1031. Elsevier, 2012.

[BBV19] Paolo Ballarini, Benoît Barbot, and Nicolas Vasselin. Performance modelling of access control mechanisms for local and vehicular wireless networks. *CoRR*, abs/1901.04285, 2019.

[BC19] Luca Bortolussi and Francesca Cairoli. Bayesian abstraction of Markov population models. In *International Conference on Quantitative Evaluation of Systems (QEST2019)*, pages 259-276. Springer, 2019.

[BCH$^{+}$07] C. Baier, L. Cloth, B. Haverkort, M. Kuntz, and M. Siegle. Model checking action- and state-labelled Markov chains. *IEEE Trans. on Software Eng.*, 33(4), 2007.

[BCMR09] Mark A. Beaumont, Jean-Marie Cornuet, Jean-Michel Marin, and Christian P. Robert. Adaptive approximate bayesian computation. *Biometrika*, 96(4):983-990, 2009.

[BCRG04] Gilles Bernot, Jean-Paul Comet, Adrien Richard, and Janine Guespin. Application of formal methods to biological regulatory networks: extending thomas' asynchronous logical approach with temporal logic. *J Theor Biol*, 229(3):339-347, Aug 2004.

[BD15] Paolo Ballarini and Marie Duflot. Applications of an expressive statistical model checking approach to the analysis of genetic circuits. *Theor. Comput. Sci.*, 599:4-33, 2015.

[BDD$^{+}$10] Paolo Ballarini, Hilal Djafri, Marie Duflot, Serge Haddad, and Nihal Pekergin. HASL: an expressive language for statistical verification of stochastic models. Technical Report TR-LACL-2010-8, LACL University of Paris-Est Créteil, 2010.

[BDD$^{+}$11a] P. Ballarini, H. Djafri, M. Duflot, S. Haddad, and N. Pekergin. COSMOS: a statistical model checker for the hybrid automata stochastic logic. In *Proceedings of the 8th International Conference on Quantitative Evaluation of Systems (QEST'11)*, pages 143-144. IEEE Computer Society Press, sep. 2011.

[BDD$^{+}$11b] Paolo Ballarini, Hilal Djafri, Marie Duflot, Serge Haddad, and Nihal Pekergin. HASL: an expressive language for statistical verification of stochastic models. In *Proc. Valuetools'2011*, pages 306-315, 2011.
[BDD<sup>+</sup>11c] Paolo Ballarini, Hilal Djafri, Marie Duflot, Serge Haddad, and Nihal Pekergin. HASL: An expressive language for statistical verification of stochastic models. In Pete Harrison Samson Lasaulce, Dieter Fiems and Luc Vandendorpe, editors, *Proceedings of the 5th International Conference on Performance Evaluation Methodologies and Tools (VALUETOOLS'11)*, pages 306–315, Cachan, France, May 2011.

[BDD<sup>+</sup>11d] Paolo Ballarini, Hilal Djafri, Marie Duflot, Serge Haddad, and Nihal Pekergin. Petri nets compositional modeling and verification of flexible manufacturing systems. In *2011 IEEE International Conference on Automation Science and Engineering*, pages 588–593, 2011.

[BDHA20] Paolo Ballarini, Davide Duma, András Horváth, and Roberto Aringhieri. Petri nets validation of markovian models of emergency department arrivals. In Ryszard Janicki, Natalia Sidorova, and Thomas Chatain, editors, *Applied and Theory of Petri Nets and Concurrency - 41st International Conference, PETRI NETS 2020, Paris, France, June 24-25, 2020, Proceedings*, volume 12152 of *Lecture Notes in Computer Science*, pages 219–238. Springer, 2020.

[BDL<sup>+</sup>12] Peter E. Bulychev, Alexandre David, Kim Guldstrand Larsen, Marius Mikucionis, Danny Bøgsted Poulsen, Axel Legay, and Zheng Wang. Uppaal-smc: Statistical model checking for priced timed automata. In Herbert Wiklicky and Mieke Massink, editors, *QAPL*, volume 85, pages 1–16, 2012.

[BDPF09] G. Balbo, M. De Pierro, and G. Franceschinis. Tagged Generalized Stochastic Petri Nets. In *Computer Performance Engineering*, volume LNCS 5652, pages 1–15, Berlin, 2009. Springer-Verlag.

[Ben21] Mahmoud Bentriou. *Statistical Inference and Verification of Chemical Reaction Networks*. PhD thesis, 12 2021.

[BG10] P. Ballarini and M.L. Guerriero. Query-based verification of qualitative trends and oscillations in biochemical systems. *Theoretical Computer Science*, 411(20):2019 – 2036, 2010.

[BGGM14] Paolo Ballarini, Emmanuelle Gallet, Pascale Le Gall, and Matthieu Manceny. Formal analysis of the wnt/β-catenin through statistical model checking. In Tiziana Margaria and Bernhard Steffen, editors, *Leveraging Applications of Formal Methods, Verification and Validation. Specialized Techniques and Applications - 6th International Symposium, ISOLA 2014, Imperial, Corfu, Greece, October 8-11, 2014, Proceedings, Part II*, volume 8803 of *Lecture Notes in Computer Science*, pages 193–207. Springer, 2014.
[BH08] Paolo Ballarini and András Horváth. Memory efficient calculation of path probabilities in large structured markov chains. In *Fifth International Conference on the Quantitative Evaluation of Systems (QEST 2008), 14-17 September 2008, Saint-Malo, France*, pages 157-166. IEEE Computer Society, 2008.

[BH09a] Paolo Ballarini and András Horváth. Compositional model checking of product-form ctmcs. *Electronic Notes in Theoretical Computer Science*, 250(1):21-37, 2009. Proceedings of the Seventh International Workshop on Automated Verification of Critical Systems (AVoCS 2007).

[BH09b] Paolo Ballarini and András Horváth. Memory efficient analysis for a class of large structured markov chains: work in progress. In Giovanni Stea, Jean Mairesse, and José Mendes, editors, *4th International Conference on Performance Evaluation Methodologies and Tools, VALUETOOLS '09, Pisa, Italy, October 20-22, 2009*, page 21. ICST/ACM, 2009.

[BH21a] Paolo Ballarini and András Horváth. Formal analysis of production line systems by probabilistic model checking tools. In *Proceedings of the 2021 IEEE Emerging Technology and Factory Automation (ETFA)*, 2021.

[BH21b] Paolo Ballarini and András Horváth. Performance analysis of production lines through statistical model checking. In Paolo Ballarini, Hind Castel, Ioannis Dimitriou, Mauro Iacono, Tuan Phung-Duc, and Joris Walraevens, editors, *Performance Engineering and Stochastic Modeling - 17th European Workshop, EPEW 2021, and 26th International Conference, ASMTA 2021, Virtual Event, December 9-10 and December 13-14, 2021, Proceedings*, volume 13104 of *Lecture Notes in Computer Science*, pages 264-281. Springer, 2021.

[BHHK00] C. Baier, B. R. Haverkort, H. Hermanns, and J.-P. Katoen. On the logical characterisation of performability properties. In *ICALP'00*, LNCS 1853, 2000.

[BHHK03a] C. Baier, B. Haverkort, H. Hermanns, and J.-P. Katoen. Model-checking algorithms for CTMCs. *IEEE Trans. on Software Eng.*, 29(6), 2003.

[BHHK03b] Christel Baier, Boudewijn Haverkort, Holger Hermanns, and Joost-Pieter Katoen. Model-checking algorithms for continuous-time markov chains. *Software Engineering, IEEE Transactions on*, 29:524-541, 07 2003.

[BHHT08] L. Bodrog, A. Heindl, G. Horváth, and M. Telek. A Markovian canonical form of second-order matrix-exponential processes. *European Journal of Operation Research*, 190:459-477, 2008.

[BHT08] L. Bodrog, A. Horváth, and M. Telek. Moment characterization of matrix exponential and Markovian arrival processes. *Annals of Operations Research*, 160:51-68, 2008.
[BKMP15] Benoît Barbot, Marta Kwiatkowska, Alexandru Mereacre, and Nicola Paolletti. Estimation and verification of hybrid heart models for personalised medical and wearable devices. In *13th International Conference on Computational Methods in Systems Biology (CMSB 2015)*, volume 9308 of LNCS, pages 3–7. Springer, 2015.

[BM06] Paolo Ballarini and Alice Miller. Model checking medium access control for sensor networks. In *Leveraging Applications of Formal Methods, Second International Symposium, ISoLA 2006, Paphos, Cyprus, 15-19 November 2006*, pages 255–262. IEEE Computer Society, 2006.

[BMM09] Paolo Ballarini, Radu Mardare, and Ivan Mura. Analysing Biochemical Oscillation through Probabilistic Model Checking. *Electronic Notes in Theoretical Computer Science*, 2009.

[BMM13] Paolo Ballarini, Lynda Mokdad, and Quentin Monnet. Modeling tools for detecting dos attacks in wsns. *Secur. Commun. Networks*, 6(4):420–436, 2013.

[BMR12] Paolo Ballarini, Jarno Mäkelä, and Andre S. Ribeiro. Expressive statistical model checking of genetic networks with delayed stochastic dynamics. In David R. Gilbert and Monika Heiner, editors, *Computational Methods in Systems Biology - 10th International Conference, CMSB 2012, London, UK, October 3-5, 2012. Proceedings*, volume 7605 of Lecture Notes in Computer Science, pages 29–48. Springer, 2012.

[BMS16] Luca Bortolussi, Dimitrios Milios, and Guido Sanguinetti. Smoothed model checking for uncertain continuous-time markov chains. *Inf. Comput.*, 247:235–253, 2016.

[BPSV05] G. Bucci, R. Piovosi, L. Sassoli, and E. Vicario. Introducing probability within state class analysis of dense-time-dependent systems. In *Second International Conference on the Quantitative Evaluation of Systems (QEST'05)*, pages 13–22, 2005.

[BPSW70] L. E. Baum, T. Petrie, G. Soules, and N. Weiss. A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. *The Annals of Mathematical Statistics*, 41(1):164–171, 1970.

[BS17] Enrico Bibbona and Roberta Sirovich. Strong approximation of density dependent Markov chains on bounded domains by jump diffusion processes. Technical report, Università di Torino, 2017.

[Buc03] P. Buchholz. *An EM-Algorithm for MAP Fitting from Real Traffic Data*, pages 218–236. Springer Berlin Heidelberg, Berlin, Heidelberg, 2003.
[BUT]

[CAH14] Marcello Colledani, Alessio Angius, and András Horváth. Lead time distribution in unreliable production lines processing perishable products. In *Proceedings of the 2014 IEEE Emerging Technology and Factory Automation (ETFA)*, pages 1–8, 2014.

[Can88] J. Canny. Some algebraic and geometric computations in pspace. In *Proc. STOC'88*, 1988.

[Car09] Luca Cardelli. *Artificial Biochemistry*, pages 429–462. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.

[CD18] José Chacón and Tarn Duong. *Multivariate Kernel Smoothing and its Applications*. 2018.

[CDM11] Taolue Chen, Marta Diciolla, Marco an Kwiatkowska, and Alexandru Mereacre. Time-bounded verification of ctmscs against real-time specifications. In *9th International Conference, FORMATS 2011*, Lecture Notes in Computer Science, pages 26–42. Springer, 2011.

[CDP$^+$17] Milan Ceska, Frits Dannenberg, Nicola Paoletti, Marta Kwiatkowska, and Lubos Brim. Precise parameter synthesis for stochastic biochemical systems. *Acta Informatica*, 54(6):589–623, 2017.

[CG08] Allan Clark and Stephen Gilmore. State-aware performance analysis with extended stochastic probes. In *Proceedings of the 5th European Performance Engineering Workshop on Computer Performance Engineering*, EPEW '08, pages 125–140, Berlin, Heidelberg, 2008. Springer-Verlag.

[CHA15] M. Colledani, A. Horváth, and Alessio Angius. Production quality performance in manufacturing systems processing deteriorating products. *Cirp Annals-manufacturing Technology*, 64:431–434, 2015.

[Che99] Song Chen. Beta kernel estimators for density functions. *Computational Statistics & Data Analysis*, 31(2):131–145, 1999.

[CHKM09] T. Chen, T. Han, J.-P. Katoen, and A. Mereacre. Quantitative model checking of CTMC against timed automata specifications. In *Proc. LICS'09*, 2009.

[CKKP05] L. Cloth, J.-P. Katoen, M. Khattri, and R. Pulungan. Model checking Markov reward models with impulse rewards. In *Proc. International Conference on Dependable Systems and Networks (DSN'05)*, pages 722–731. IEEE Computer Society Press, 2005.
[COS] Cosmos home page. https://cosmos.lacl.fr/.

[CP34] C. J. Clopper and E. S. Pearson. The use of confidence or fiducial limits illustrated in the case of the binomial. *Biometrika*, 26(26):404-413, 1934.

[CR65] Y. S. Chow and Herbert Robbins. On the asymptotic theory of fixed-width sequential confidence intervals for the mean. *Annals of Mathematical Statistics*, 36(2):457-462, 1965.

[CRCD+04] Nathalie Chabrier-Rivier, Marc Chiaverini, Vincent Danos, François Fages, and Vincent Schächter. Modeling and querying biomolecular interaction networks. *Theoretical Computer Science*, 325(1):25-44, 2004. Computational Systems Biology.

[CT12] M. Colledani and T. Tolio. Integrated quality, production logistics and maintenance analysis of multi-stage asynchronous manufacturing systems with degrading machines. *CIRP Annals-Manufacturing Technology*, 61/1:455-458, 2012.

[DHK04] Nicholas J. Dingle, Peter G. Harrison, and William J. Knottenbelt. Uniformisation and Hypergraph Partitioning for the Distributed Computation of Response Time Densities in Very Large Markov Models. *Journal of Parallel and Distributed Computing*, 64(8):309-920, August 2004.

[DHS07] Susanna Donatelli, Serge Haddad, and Jeremy Sproston. Csl<sup>TA</sup>: an expressive logic for continuous-time markov chains. In *Fourth International Conference on the Quantitative Evaluation of Systems (QEST 2007)*, pages 31-40, 2007.

[DHS09] S. Donatelli, S. Haddad, and J. Sproston. Model checking timed and stochastic properties with CSL<sup>TA</sup>. *IEEE Trans. on Software Eng.*, 35, 2009.

[DJKV17] Christian Dehnert, Sebastian Junges, Joost-Pieter Katoen, and Matthias Volk. A storm is coming: A modern probabilistic model checker. In *Proc. 29th International Conference on Computer Aided Verification (CAV'17)*, 2017.

[DM10] Alexandre Donzé and Oded Maler. Robust satisfaction of temporal logic over real-valued signals. In *Proceedings of the 8th International Conference on Formal Modeling and Analysis of Timed Systems, FORMATS'10*, pages 92-106, Berlin, Heidelberg, 2010. Springer-Verlag.

[DMDJ06] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential monte carlo samplers. *Journal of the Royal Statistical Society B*, 68(3):411-436, 2006.

[EC80] E. A. Emerson and E. M. Clarke. Characterizing correctness properties of parallel programs using fixpoints. In *Proc. of the 7th Coll. on Automata, Languages and Programming*, LNCS 85, 1980.
[EH86] E. Allen Emerson and Joseph Y. Halpern. "sometimes" and "not never" revisited: on branching versus linear time temporal logic. *J. ACM*, 33(1):151-178, 1986.

[EL00] Michael B. Elowitz and Stanislas Leibier. A synthetic oscillatory network of transcriptional regulators. *Nature*, 2000.

[Gil77] D.T. Gillespie. Exact stochastic simulation of coupled chemical reactions. *Journal of Physical Chemistry*, 81(25):2340-2361, 1977.

[Gil01] Daniel T. Gillespie. Approximate accelerated stochastic simulation of chemically reacting systems. *The Journal of Chemical Physics*, 115(4):1716-1733, 2001.

[Gly83] P. W. Glynn. On the role of generalized semi-Markov processes in simulation output analysis. In *Proc. Conf. Winter simulation*, 1983.

[GMGB14] Emmanuelle Gallet, Matthieu Manceny, Pascale Le Gall, and Paolo Ballarini. An LTL model checking approach for biological parameter inference. In Stephan Merz and Jun Pang, editors, *Formal Methods and Software Engineering - 16th International Conference on Formal Engineering Methods, ICFEM 2014, Luxembourg, Luxembourg, November 3-5, 2014. Proceedings*, volume 8829 of *Lecture Notes in Computer Science*, pages 155-170. Springer, 2014.

[GMGB15] Emmanuelle Gallet, Matthieu Manceny, Pascale Le Gall, and Paolo Ballarini. Étude de réseaux de thomas par validation de propriétés LTL pour pseudomonas aeruginosa. *Tech. Sci. Informatiques*, 34(5):575-600, 2015.

[Gol02] Albert Goldbeter. Computational approaches to cellular rhythms. *Nature*, 420:238-245, 2002.

[HHWT97] Thomas A. Henzinger, Pei-Hsin Ho, and Howard Wong-Toi. Hytech: A model checker for hybrid systems. In Orna Grumberg, editor, *Computer Aided Verification*, pages 460-463, Berlin, Heidelberg, 1997. Springer Berlin Heidelberg.

[Hil05] J. Hillston. Process algebras for quantitative analysis. In *Proceedings of the 20th Annual IEEE Symposium on Logic in Computer Science*, pages 239-248, Washington, DC, USA, 2005. IEEE Computer Society.

[HJ90] Hans Hansson and Bengt Jonsson. A logic for reasoning about time and reability. Technical report, 1990.

[HJB$^{+}$10] R. He, P. Jennings, S. Basu, A. P. Ghosh, and H. Wu. A bounded statistical approach for model checking of unbounded until properties. In *Proc. ASE'10*, 2010.
[HJK+22] Christian Hensel, Sebastian Junges, Joost-Pieter Katoen, Tim Quatmann, and Matthias Volk. The probabilistic model checker storm. *International Journal on Software Tools for Technology Transfer*, 24(4):589–610, 2022.

[HKM08] T. Han, J. Katoen, and A. Mereacre. Approximate parameter synthesis for probabilistic time-bounded reachability. In *2008 Real-Time Systems Symposium*, pages 173–182, 2008.

[HKPV98] Thomas A. Henzinger, Peter W. Kopke, Anuj Puri, and Pravin Varaiya. What's decidable about hybrid automata? *Journal of Computer and System Sciences*, 57(1):94–124, 1998.

[HLMP04] Thomas Hérault, Richard Lassaigne, Frédéric Magniette, and Sylvain Peyronnet. Approximate probabilistic model checking. In Bernhard Steffen and Giorgio Levi, editors, *Verification, Model Checking, and Abstract Interpretation*, pages 73–84, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg.

[HLP06] T. Herault, R. Lassaigne, and S. Peyronnet. APMC 3.0: Approximate verification of discrete and continuous time Markov chains. In *Proc. QEST'06*, 2006.

[Hoe63] Wassily Hoeffding. Probability inequalities for sums of bounded random variables. *Journal of the American Statistical Association*, 58(301):pp. 13–30, 1963.

[Hor13] G. Horváth. Matching marginal moments and lag autocorrelations with maps. In *Proceedings of the 7th International Conference on Performance Evaluation Methodologies and Tools*, ValueTools '13, pages 59–68, ICST, Brussels, Belgium, Belgium, 2013. ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering).

[HPRV11] András Horváth, Marco Paolieri, Lorenzo Ridi, and Enrico Vicario. Probabilistic model checking of non-markovian models with concurrent generally distributed timers. In *Eighth International Conference on Quantitative Evaluation of Systems, QEST 2011, Aachen, Germany, 5-8 September, 2011*, pages 131–140. IEEE Computer Society, 2011.

[HRS13] M Heiner, C Rohr, and M Schwarick. MARCIE - Model checking And Reachability analysis done effiCiently. In JM Colom and J Desel, editors, *Proc. PETRI NETS 2013*, volume 7927 of LNCS, pages 389–399. Springer, June 2013.

[HT07] A. Horváth and M. Telek. Matching more than three moments with acyclic phase type distributions. *Stochastic Models*, 23(2):167–194, 2007.
[JLS12] Cyrille Jégourel, Axel Legay, and Sean Sedwards. A platform for high performance statistical model checking - plasma. In Cormac Flanagan and Barbara König, editors, *TACAS*, volume 7214, pages 498–503. Springer, 2012.

[JSD18] Cyrille Jegourel, Jun Sun, and Jin Song Dong. On the sequential massart algorithm for statistical model checking. In Tiziana Margaria and Bernhard Steffen, editors, *Leveraging Applications of Formal Methods, Verification and Validation. Verification*, pages 287–304, Cham, 2018. Springer International Publishing.

[KBVC16] Konstantinos Koutroumpas, Paolo Ballarini, Irene Votsi, and Paul Henry Cournède. Bayesian parameter estimation for the Wnt pathway: An infinite mixture models approach. volume 32, pages 781–789, 2016.

[KCAB09] Zohra Khalis, Jean-Paul Comet, Richard Adrien, and Gilles Bernot. The smbionet method for discovering models of gene regulatory networks. *Genes, Genomes and Genomics*, 3, 01 2009.

[KM27] W. O. Kermack and A. G. McKendrick. A Contribution to the Mathematical Theory of Epidemics. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 115(772):700–721, 1927.

[KNP02] Marta Kwiatkowska, Gethin Norman, and David Parker. Prism: Probabilistic symbolic model checker. pages 200–204. Springer, 2002.

[KNP07a] M. Kwiatkowska, G. Norman, and D. Parker. Stochastic model checking. In *Formal Methods for the Design of Computer, Communication and Software Systems: Performance Evaluation*, volume 4486 of LNCS, pages 220–270. Springer, 2007.

[KNP07b] Marta Z. Kwiatkowska, Gethin Norman, and David Parker. Stochastic model checking. In Marco Bernardo and Jane Hillston, editors, *SFM*, volume 4486, pages 220–270. Springer, 2007.

[KNP11] M. Kwiatkowska, G. Norman, and D. Parker. PRISM 4.0: Verification of probabilistic real-time systems. In G. Gopalakrishnan and S. Qadeer, editors, *Proc. 23rd International Conference on Computer Aided Verification (CAV'11)*, volume 6806 of LNCS, pages 585–591. Springer, 2011.

[KNS02] Marta Kwiatkowska, Gethin Norman, and Jeremy Sproston. Probabilistic model checking of the ieee 802.11 wireless local area network protocol. In Holger Hermanns and Roberto Segala, editors, *Process Algebra and Probabilistic Methods: Performance Modeling and Verification*, pages 169–187, Berlin, Heidelberg, 2002. Springer Berlin Heidelberg.
[Knu97] Donald E. Knuth. *The Art of Computer Programming, Volume 2 (3rd Ed.): Seminumerical Algorithms*. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1997.

[Kul95] V.G. Kulkarni. *Modeling and Analysis of Stochastic Systems*. Chapman & Hall, London, UK, 1995.

[Kul16] V.G. Kulkarni. *Modeling and Analysis of Stochastic Systems, Third Edition*. Chapman & Hall/CRC Texts in Statistical Science. CRC Press, 2016.

[Kur76] Thomas G Kurtz. Limit theorems and diffusion approximations for density dependent Markov chains. In *Stochastic Systems: Modeling, Identification and Optimization, I*, pages 67–78. Springer, 1976.

[KZH$^{+}$09] Joost-Pieter Katoen, Ivan S. Zapreev, Ernst Moritz Hahn, Holger Hermanns, and David N. Jansen. The ins and outs of the probabilistic model checker MRMC. *QEST*, pages 167–176, 2009.

[LBGY16] Christopher Lester, Ruth E Baker, Michael B Giles, and Christian A Yates. Extending the multi-level method for the simulation of stochastic biological systems. *Bulletin of mathematical biology*, 78(8):1640–1677, 2016.

[LKKP05] L.Cloth, J.-P. Katoen, M. Khattri, and R. Pulungan. Model checking Markov reward models with impulse rewards. In *Proc. DSN’05*, 2005.

[Lot09] Alfred J. Lotka. Contribution to the theory of periodic reactions, January 1909.

[LST16] Axel Legay, Sean Sedwards, and Louis-Marie Traonouez. Plasma lab: A modular statistical model checking platform. In *ISoLA (1)*, volume 9952 of *Lecture Notes in Computer Science*, pages 77–93, 2016.

[MMB$^{+}$17] Quentin Monnet, Lynda Mokdad, Paolo Ballarini, Youcef Hammal, and Jalel Ben-Othman. Dos detection in wsns: Energy-efficient methods for selecting monitoring nodes. *Concurr. Comput. Pract. Exp.*, 29(23), 2017.

[MN04] Oded Maler and Dejan Nickovic. Monitoring temporal properties of continuous signals. In Yassine Lakhnech and Sergio Yovine, editors, *Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems*, pages 152–166, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg.

[MPRR12] Jean-Michel Marin, Pierre Pudlo, Christian P. Robert, and Robin J. Ryder. Approximate bayesian computational methods. *Statistics and Computing*, 22(6):1167–1180, Nov 2012.
[Neu75] M.F. Neuts. Probability distributions of phase type. In *Liber Amicorum Prof. Emeritus H. Florin*, pages 173–206. University of Louvain, 1975.

[Neu79] M.F. Neuts. A versatile markovian point process. *Journal of Applied Probability*, 16:764–779, 1979.

[OD09] H. Okamura and T. Dohi. Faster maximum likelihood estimation algorithms for markovian arrival processes. In *Quantitative Evaluation of Systems, 2009. QEST'09. Sixth International Conference on the*, pages 73–82. IEEE, 2009.

[OS99] W. Douglas Obal, II and William H. Sanders. State-space support for path-based reward variables. *Perform. Eval.*, 35:233–251, May 1999.

[Pnu77] Amir Pnueli. The temporal logic of programs. In *Proceedings of the 18th Annual Symposium on Foundations of Computer Science*, SFCS '77, pages 46–57, Washington, DC, USA, 1977. IEEE Computer Society.

[PRI] Prism home page. http://www.prismmodelchecker.org.

[QS82] J. Queille and J. Sifakis. Specification and verification of concurrent systems in cesar. pages 337–351. 1982.

[Rab90] L. R. Rabiner. Readings in speech recognition. chapter A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition, pages 267–296. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1990.

[RP09] Diana El Rabih and Nihal Pekergin. Statistical model checking for steady state dependability verification. In *2009 Second International Conference on Dependability*, pages 166–169, 2009.

[Ryd96] T. Ryden. An em algorithm for estimation in markov-modulated poisson processes. *Computational Statistics & Data Analysis*, 21(4):431 – 447, 1996.

[RZK06] A. Ribeiro, R. Zhu, and S. A. Kauffman. A general modeling strategy for gene regulatory networks with stochastic dynamics. *Journal of computational biology: a journal of computational molecular cell biology*, 13(9):1630–1639, November 2006.

[SFB18] Scott A Sisson, Yanan Fan, and Mark Beaumont. *Handbook of approximate Bayesian computation*. Chapman and Hall/CRC, 2018.

[Sil86] B W Silverman. *Density Estimation for Statistics and Data Analysis*. Chapman & Hall, London, 1986.

[Spi13] D. Spieler. Characterizing oscillatory and noisy periodic behavior in markov population models. In *Proc. QEST'13*, 2013.
[SVA05a] K. Sen, M. Viswanathan, and G. Agha. VESTA: A statistical model-checker and analyzer for probabilistic systems. In Proc. *QEST'05*, 2005.

[SVA05b] Koushik Sen, Mahesh Viswanathan, and Gul Agha. On statistical model checking of stochastic systems. In *CAV*, volume 3576 of *Lecture Notes in Computer Science*, pages 266–280. Springer, 2005.

[TH07] M. Telek and G. Horváth. A minimal representation of Markov arrival processes and a moments matching method. *Performance Evaluation*, 64(9-12):1153–1168, Aug. 2007.

[VKBL02] J. Vilar, H.-Y. Kueh, N. Barkai, and S. Leibler. Mechanisms of noise-resistance in genetic oscillators. *Proc. National Academy of Sciences of the United States of America*, 99(9):5988–5992, 2002.

[Vol28] Vito Volterra. Variations and Fluctuations of the Number of Individuals in Animal Species living together. *ICES Journal of Marine Science*, 3(1):3–51, 04 1928.

[Wal45] A. Wald. Sequential tests of statistical hypotheses. *The Annals of Mathematical Statistics*, 16(2):117–186, 06 1945.

[You05] Håkan LS Younes. Ymer: A statistical model checker. In *Computer Aided Verification*, pages 429–433. Springer, 2005.

[YS06] H.L.S. Younes and R.G. Simmons. Statistical probabilistic model checking with a focus on time-bounded properties. *Inf. Comput.*, 204(9), 2006.

[ZC13] André Zuliani, Paolo an Platzer and Edmund M. Clarke. Bayesian statistical model checking with application to stateflow/simulink verification. *Formal Methods in System Design*, 43(2):338–367, 2013.