Christophe

Xiangtuo

Erwan, Sylvie), ceux arrivés après (Walid, Théo, Gurvan Alexandre

Agathe Laura

Tout d'abord, je tenais à remercier mes directeurs de thèse Paul-Henry Cournède et Paolo Ballarini pour l'opportunité de réaliser cette thèse stimulante par ses différents aspects complèxes et variés. Merci Paul-Henry de m'avoir accordé ta confiance dés le début de notre rencontre, d'avoir su être à l'écoute durant toutes ces années pour les concrétiser en l'obtention de mon doctorat, en simultanéité avec ta promotion de Directeur de la Recherche de CentraleSupélec (encore félicitations). Merci aux rapporteur.trice.s de thèse Adeline Leclercq-Samson et Blaise Genest d'avoir accepté de faire parti du jury, tant pour leurs retours pertinents sur le manuscrit que les questions soulevées impregnées d'intérêt durant la soutenance. Merci à Benoît Barbot qui a rapidement accepté d'examiner le manuscrit après une indisponibilité d'un des membres du jury. Ces années auraient été bien plus tristes sans l'ambiance et la bienveillance des membres du laboratoire MICS. J'ai été plus que chanceux d'avoir partagé mon bureau avec ces deux personnes brillantes et drôles que sont Brice et Antonin. Bien entendu, je remercie la team Kwak pour les afterworks de qualité. Mathilde, on a commencé ensemble, désolé d'avoir pris un peu de retard. Elvire, bon courage pour la dernière ligne droite. Merci à l'équipe du mésocentre (Laurent, Rémi, Guillaume) qui a permis de rendre le travail de cette thèse bien plus agréable. Merci à Sylvain pour tes meilleures recommandations de festivals, Gautier pour tes tentatives journalières manquées de me faire rendre le manuscrit à temps

Résumé

Les réseaux de réactions chimiques (CRN) constituent un formalisme utilisé pour modéliser des processus biologiques. Quand la population est de taille modérée et supposée bien mélangée, le processus stochastique sous-jacent pour décrire ses dynamiques est une chaîne de Markov en temps continu (CTMC). Ce processus est dit sans mémoire: l'état futur du système ne dépend que de l'état courant.

L'inférence statistique de ce type de CTMC est complexe: le calcul de la vraisemblance est en général difficile à résoudre. Les méthodes ABC (Approximate Bayesian Computation) forment une classe de méthodes bayésiennes sans calcul de vraisemblance qui permettent d'approcher la distribution postérieure avec des simulations de Monte-Carlo.

La vérification de modèles, qui fut à l'origine développée pour garantir la fiabilité de systèmes et logiciels informatiques, se penche de plus en plus sur la biologie des systèmes. En effet, il y a un réel besoin de comprendre les interactions complexes entre molécules dans les systèmes biologiques. Malheureusement, l'espace d'états d'un CTMC défini par un CRN explose généralement, voir est infini. Pour palier à cela, des méthodes de vérification statistiques ont été développées. Le principe est de simuler un certain nombre de fois le modèle et de calculer le ratio des simulations qui ont vérifié une propriété. Récemment, une logique temporelle appelée HASL a été introduite pour la vérification statistique de modèles: elle adopte intrinsèquement le point de vue statistique de la vérification.

Dans cette thèse, nous nous intéressons à l'inférence statistique et la vérification statistique de chaines de Markov en temps continu définies par un modèle de réseaux de réactions chimiques. Notre contribution tient principalement dans la formulation d'un algorithme ABC combiné avec le formalisme HASL appelé automaton-ABC. xiv

List of Figures

2.2

Graph of the SIR CTMC with the initial state s 0 = (95, 5, 0). . . . . .

2.3

Effect of the population size on the simulations of the SIR model (only the number of infected people is plotted) with θ = (0.12/N pop , 0.05) and τ = 5.0. On the left plot, s 0 = (95, 5, 0) (N pop = 100). On the right plot, s 0 = (9500, 500, 0) (N pop = 10000). Red paths are simulated from the SSA. Green paths are simulated from the Tau-leaping algorithm.

5 paths for each algorithm. . . . . . . . . . . . . . . . . . . . . . . . .

2.4

Effect of the tau-leap time step on the simulations of the SIR model (only the number of infected people is plotted) with θ = (0.12/N pop , 0.05) and s 0 = (95, 5, 0) (N pop = 100). In the left plot, τ = 0.5. In the right plot, τ = 5.0. Red paths are simulated from the SSA. Green paths are simulated from the Tau-leaping algorithm. 5 paths for each algorithm. 5.12 5 simulated trajectories of the species P 1 . β varies in {0.5, 1.0, 2.0, 4.0}.

5.13 5 simulated trajectories of the species P 1 . n varies in {0.5, 1.0, 2.0, 5.0}.

5.14 5 simulated trajectories of the species P 1 . α 0 varies in {0.0, 0.01, 0.1, 1.0}.103

Correlation plot of automaton-ABC posterior with

6.13

Results for the SIR model with φ = G [0,100] (I > 0) ∧ F [100,120] (I = 0). Markov Chains [[p, q]] The set of integers {p, . . . , q -1, q} The kinetic rate of the reaction channel R j ν j

The stoichiometric vector of the reaction channel R j

Statistical observations

i.i.d independent and identically distributed 1 : N The set {1, . . . , N } X (1) , . . . , X (N ) N i.i.d random observations

x (1) , . . . , x (N ) A realisation of N i.i. and the media has seized upon two types of questions only scientists were interested in so far:

1. How can one infer the spread of the disease based on previously infected people data?

2. How can we specify properties such as in which region vigilance should be reinforced, the epidemic is stable or over?

The first problem is related to statistical inference [START_REF] Casella | Statistical Inference[END_REF].

Based on probability theory, statistical inference constructs estimators of properties of a general population from samples, quantifies the uncertainty of such estimation and predicts possible future observations based on available data.

The second problem is related to model checking [START_REF] Baier | Principles of model checking[END_REF].

Based on theoretical computer science, this field groups techniques to check if a model fulfils some specification φ expressed by formal logic. It provides methods that explore the state-space model efficiently in order to assess whether the model verifies the specification φ or not. Initially, it was developed to assess hardware and software systems' reliability, but there is a growing interest in model checking of Systems Biology [START_REF] Kwiatkowska | Using Probabilistic Model Checking in Systems Biology[END_REF].

In this thesis, we focus on statistical inference and verification of stochastic chemical kinetics. Chemical Reaction Network is a formalism that provides a way to model the interactions between entities. A reaction is an event that can occur in the system, 1.1. Context which modifies each entity's quantities. Even if this formalism is originally used to describe chemical systems [START_REF] Feinberg | Foundations of Chemical Reaction Network Theory[END_REF], it is also used in gene regulatory networks [START_REF] Karlebach | Modelling and analysis of gene regulatory networks[END_REF] or systems of interactions between humans. An illustrative example that is a guideline model in the thesis is the SIR epidemiology model:

Inf ection : S + I k i -→ 2I Recovery : I kr -→ R (1.1)
This Chemical Reaction Network has three types of entities (living beings in this context): S is the population that can be infected, I is the infected people, whereas R is the recovered/removed population. The system contains two reactions: either an infection or an immunisation/death can occur.

Chemical Reaction Networks are algebraic models, but they induce mathematical models. Traditionally, their dynamics are described by a system of Ordinary Differential Equations (ODE). Such kind of dynamics is continuous and deterministic. If the deterministic model is relevant under some assumptions, such as a high enough species density, they are not adapted to populations of small sizes. With relaxed assumptions, the stochastic dynamics are described by a Continuous-Time Markov Chain (CTMC). The main property of such a stochastic process is the memoryless property (also known as Markov property): the system's future state only depends on the current state. For example, genetic networks have proven their stochastic nature [START_REF] Thattai | Intrinsic noise in gene regulatory networks[END_REF].

Statistical inference of CTMCs has been widely studied in the literature, particularly Bayesian statistics [START_REF] Craciun | Statistical model for biochemical network inference[END_REF][START_REF] Schnoerr | Approximation and inference methods for stochastic biochemical kinetics -A tutorial review[END_REF][START_REF] Loskot | Comprehensive review of models and methods for inferences in bio-chemical reaction networks[END_REF]. A general issue raised by such a complex stochastic model is the evaluation of the likelihood: it is intractable in most cases. The Bayesian framework allows the use of preexisting expertise on models via the prior distribution and quantifies the uncertainty of parameter inference with the posterior distribution. Also, it offers groups of likelihood-free methods, which allows approximating the posterior distribution. Approximate Bayesian Computation, a specific family of likelihood-free methods, is specifically studied in this thesis as our models have intractable likelihoods. It relies on Monte Carlo simulation methods:

the general idea is to draw a parameter θ, simulate the model y ∼ p(.|θ), and keep the sample if y is close enough to observations y exp . Since its first appearance in the literature of population genetics [START_REF] Tavaré | Inferring coalescence times from DNA sequence data[END_REF][START_REF] Pritchard | Population Growth of Human Y Chromosomes: A Study of Y Chromosome Microsatellites[END_REF] ABC methods have received much attention from the statistician community [START_REF] Marin | Approximate Bayesian Computational methods[END_REF]) [START_REF] Sisson | Overview of Approximate Bayesian Computation[END_REF]. They have proven efficiency in the case of CTMCs [START_REF] Warne | Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art[END_REF][START_REF] Alharbi | Bayesian Inference for Continuous Time Markov Chains[END_REF].

Chapter 1. Introduction 3 Verification of a CTMC M has been studied for about two decades since the appearance of CSL temporal logic [START_REF] Aziz | Verifying Continuous Time Markov Chains[END_REF]. As CTMCs are stochastic, verification of a CTMC considers the probability of verifying a specification φ. This kind of methods is divided into two approaches. First, the qualitative approach assesses whether the probability is above (or below) a threshold. Second, the quantitative approach is related to the estimation of such probabilities. Numerical methods have been developed for these purposes [START_REF] Baier | Model-Checking Algorithms for Continuous-Time Markov Chains[END_REF]. To face the problem of state space explosion, statistical methods for verification have been considered since 2002 [START_REF] Legay | Statistical model checking: An overview[END_REF][START_REF] Legay | Statistical Model Checking[END_REF]. This was enhanced by the growing interest of Systems Biology field in verification methods: the state-space of CTMCs defined by a Chemical Reaction Network quickly explodes or is infinite.

Initially, these methods were developed for a single CTMC, but now verification methods tackle model checking with CTMCs indexed by parameters (parametric CTMC).

The qualitative approach for parametric CTMC is called parameter synthesis. An approximate numerical method has been first proposed [START_REF] Han | Approximate parameter synthesis for probabilistic time-bounded reachability[END_REF], and also statistical methods [START_REF] Bortolussi | Bayesian statistical parameter synthesis for linear temporal properties of stochastic models[END_REF]. For the quantitative approach, statistical methods have been developed [START_REF] Bortolussi | Smoothed model checking for uncertain Continuous-Time Markov Chains[END_REF].

Outline

This thesis focuses on statistical inference and verification of parametric Continuous-Time Markov Chains defined by Chemical Reaction Networks. Our main methodological contribution relies on the formulation of a general algorithm called automaton-ABC, which uses HASL formalism within the ABC method. With this high-level algorithm, we address several tasks of statistical inference and statistical verification.

The considered tools and problems are both issued from model checking and statistical inference, which positions this thesis at the frontier between these two fields.

Efforts have been made to be intelligible for scientists of both fields.

The first three chapters are related to the literature.

• Chapter 2 presents the basics of Markov Chains theory. The main properties of Continuous-Time Markov Chains (CTMC) are discussed. We define the probability measure over the set of càdlàg functions based on a CTMC. Then, Chemical Reaction Network formalism is defined. We see how a CTMC describes the stochastic dynamics, we discuss the other possible equations of the The last two chapters expose contributions.

• Chapter 5 details the statistical inference framework of CTMCs. We present the automaton-ABC algorithm, which combines the synchronised simulation offered by HASL formalism and ABC. We present two applications based on automaton-ABC: detection of oscillatory trends in parametric CTMC and faster implementation of classical ABC inference.

• Chapter 6 details a new statistical method based on automaton-ABC for the quantitative approach of model checking parametric CTMC in time-bounded reachability problems.

The new methods lead to a package in the cross-platform Julia language [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF]. The implementation is detailed in Appendix A. Efforts have been made to make it easy to use, efficient and tested. The Scripts List showed above in the thesis groups the scripts that mainly produce the results based on the package. In Section 2.1.1, we recall definitions and basic properties of Markov Chains in a classical way for stochastic processes. Section 2.1.2 gives a Markov Chain description frequently used in model checking and based on the transition rate matrix. We describe the set of paths of a CTMC. In Section 2.2, we construct a probability measure over the set of paths of a CTMC. This section is a requisite for the definition of verification tasks of CTMC. In the last Section 2.3, we describe Chemical Reaction Networks formalism for the description of chemical kinetics. We show that the underlying probabilistic model that rules the stochastic dynamics is a CTMC. We 2.1. Different perspectives on Markov Chains also detail the equations that rule the dynamics, under new assumptions or not, and expose simulation algorithms of trajectories.

The notations related to this chapter are described in the Notations part at the beginning of the thesis. We assume the reader is familiar with classical probability and measure theory notions such as random variable, σ-algebra and probability measure.

For further exploration of the theory of Markov Chains, we refer the reader to the books [START_REF] Bladt | Matrix-Exponential Distributions in Applied Probability[END_REF], [START_REF] Stroock | An Introduction to Markov Processes[END_REF], [START_REF] Kulkarni | Modeling and Analysis of Stochastic Systems[END_REF].

In the whole chapter, Ω is the sample space with ω ∈ Ω an elementary event, F a σ-field of Ω and (Ω, F, P) a probability space. F X is the σ-algebra generated by

a set X. S = {(s (i) ) i∈I } ⊆ N d with I ⊆ N

Different perspectives on Markov Chains

This section aims to define discrete space Markov Chains in discrete time and continuous time and detail their main properties. Proof and details of these results can be found in [START_REF] Bladt | Matrix-Exponential Distributions in Applied Probability[END_REF].

Markov Chains as a stochastic process

Discrete-Time Markov Chain (DTMC)

A Discrete-Time Markov Chain is a stochastic process indexed by an integer (T = N) that fulfils the Markov property. This property states that the current state of the chain only depends on the last state of the chain.

Definition 2.1.1 (Discrete-Time Markov Chain)

(S n ) n∈N is a time-homogeneous Discrete-Time Markov Chain with state space S ⊆ N d if 1. ∀n ∈ N, ∀s, s n , . . . , s 0 ∈ S, P(S n+1 = s|S n = s n , . . . , S 0 = s 0 ) = P(S n+1 = s|S n = s n ) (Markov property) 2. ∀n ∈ N, ∀s, s ′ ∈ S, P(S n+1 = s ′ | S n = s) = P(S 1 = s ′ |S 0 = s) (Time- homogeneity)
Chapter 2. Markov Chains and Chemical Reaction Networks 7

In the following, we introduce some classical notations about Markov Chains.

They are involved in helpful results such as ergodic theorems.

Transition probabilities

Definition 2.1.2 (Transition matrix for a DTMC)

For a time-homogeneous DTMC (S n ) n∈N with state space S = {(s (i) ) i } ⊆ N d , the matrix

P = (p ij ) ij = (P(S 1 = s (j) |S 0 = s (i) )) s (i) ,s (j) ∈S
is called the transition matrix of the DTMC. P (s (i) , s (j) ) = p ij is called a transition probability.

p ij is the probability that, being at the state s (i) , the chain transitions to the state s (j) . One can generalise this definition with n-step transitions.

Definition 2.1.3 (n-step transition probability)

Let (S n ) n∈N be a time-homogeneous DTMC. For n ∈ N * , the probability p (n) ij = P(S n = s (j) |S 0 = s (i) )
is called a n-step transition probability. We denote Let (S n ) n∈N be a time-homogeneous DTMC with associated transition matrices (P (n) ) n∈N

P (n) = (p (n) ij ) ij .

Any transition probability p

Then:

∀n, m ∈ N, P (n+m) = P (n) P (m) .

Ergodic chain A stationary distribution is a distribution invariant by the transition probability matrix (see Definition 2.1.10 further). The existence of such distributions is an important result of Markov Chains. Under adequate hypotheses, the matrix P (n) converges to the unique stationary distribution of the Markov Chain as n grows to infinity. These results are known as ergodic theorems. They are especially useful in advanced Bayesian statistical methods, as we will see in the next chapter.

We proceed by describing these hypotheses related to the states of the chain. We say that a chain verifies a property if all its states verify this property.

Different perspectives on Markov Chains Definition 2.1.4 (Initial passing time)

The random time of first entrance in a state s (i) , also called initial passing time or first hit time, is defined as:

τ s (i) = inf{t ∈ N, t > 1/S t = s (i) } Definition 2.1.5 (Recurrent state) A state s (i) is recurrent if +∞ n=1 p (n) ii = +∞. Otherwise it is called transient.
This means that as n tends to infinity, the state will be visited by the Markov Chain regularly. Indeed +∞ n=1 p

(n)
ii is the expectation of the random variable N s (i) that counts the number of times the chain sojourns in s (i) when starting from s (i) , i.e.

N s (i) = +∞ n=1 1(S n = s (i) ).
Definition 2.1.6 (Positive recurrent state)

A recurrent state s (i) is positive recurrent if E[τ s (i) | S 0 = s (i) ] < +∞. Otherwise it is called null recurrent.
A Markov Chain is positive recurrent if all its states are positive recurrent.

Definition 2.1.7 (Communication of states -Irreducibility) i) and s (j) communicate (s (i) ↔ s (j) ) if s (i) → s (j) and s (j) → s (i) .

A state s (i) leads to s (j) (s (i) → s (j) ) if ∃n ∈ N, p (n) ij > 0. States s (
A Markov Chain is irreducible if all its states communicate. Definition 2.1.8 (Period of a state and aperiodicity)

The period of a state s (i) is defined as:

gcd{n ≥ 1/p (n) ii > 0}
where gcd is the greatest common divisor. If the period is 1, s (i) is called aperiodic.

A Markov Chain is aperiodic if all its states are aperiodic.

We obtain the notion of ergodicity by gathering all these properties. Definition 2.1.9 (Ergodic Markov Chain)

A Markov Chain (S n ) n∈N is ergodic if it is positive reccurent, irreducible and aperiodic.

Ergodic Markov Chains have useful properties, which are mainly embodied by the ergodic theorem.

Definition 2.1.10 (Stationary distribution)

A stationary distribution (π i ) i∈N/s (i) ∈S of a Markov Chain with transition matrix P is a non-zero vector where for all i ∈ N corresponding to a state s

(i) ∈ S, 0 ≤ π i < 1, π = πP and i∈N/s (i) ∈S π i = 1. Proposition 2.1.1 (Stationary distribution of a Markov Chain)
An ergodic Markov Chain admits a unique stationary distribution π.

In fact, it is sufficient for the Markov Chain to be irreducible and positive recurrent for Proposition 2.1.1 to be true. We have stated the required definitions and properties for the ergodic theorem.

|p (n) ij -π j | → n→+∞ 0.
This result is generalisable with continuous state spaces. They are especially useful in Markov Chain Monte Carlo statistical methods, as we will see in Chapter 3.

Continuous-Time Markov Chain

We now extend the definitions described in the previous section to the case of a continuous-time set T = R ≥0 (state space S is still discrete).

Definition 2.1.11 (Continuous-Time Markov Chain (CTMC))

(S t ) t∈R ≥0 is a time-homogeneous CTMC if it satisfies the following properties: 1. The trajectories (t → S t (ω)) ω∈Ω are right-continuous 2. ∀t, s ∈ T, P(S t |S v , v ∈ [0, s]) = P(S t |S s ) (Markov property) 3. ∀t, v ∈ R ≥0 , t > v, P(S t |S v ) = P(S t-v |S 0 ) (Time-homogeneity).
Remark 2.1.1 (Time-homogeneity of a Markov Chain)

In our applications, unless otherwise stated, Markov Chains will be considered timehomogeneous.

Different perspectives on Markov Chains

We define T s (s ∈ S) the random sojourn time in state s, i.e. the time before a state transition occurs while in s:

P(T s ≥ t) = P({ω ∈ Ω, ∀t ′ ∈ [0, t], S t ′ (ω) = s}) Proposition 2.1.2 (Exponential times of a CTMC.)
The random variable T s verifies the memoryless property, i.e. it follows an exponential law.

As T s follows an exponential law, we denote E(s) its parameter, which is called

total exit rate. Then, T s ∼ Exp(E(s)) and E(T s ) = 1 E(s)
. The event {S Ts = s ′ } is defined as the set of elementary events ω so that the trajectory S (ω) reaches s ′ after a jump from s.

Definition 2.1.12 (Transition probability for a CTMC)

Let (S t ) t∈R ≥0 be a CTMC. ∀s, s ′ ∈ S, s ̸ = s ′ , P (s, s ′ ) = P(S Ts = s ′ |S 0 = s) is called a transition probability.
Then, we can prove the following properties based on the first-order Taylor's expansion of exponential and on the Markov property:

∀t ∈ T, ∀s, s ′ ∈ S, s ̸ = s ′ ,P(S t+h = s ′ |S t = s) = P (s, s ′ )E(s)h + o(h) ∀t ∈ T, ∀s ∈ S,P(S t+h = s|S t = s) = 1 - s ′ ̸ =s P (s, s ′ )E(s)h + o(h).
(2.1)

If we are in state s at t, the probability that the Markov process jumps to s ′ after a time h only depends on h and s ′ , but does not depend on t. Let {S Ts 0 = s ′ } the event where trajectories go to s ′ after a jump from s 0 . We obtain the following result:

Proposition 2.1.3
For s ∈ S and I ⊆ R ≥0 , the events {S Ts = s ′ } and {T s ∈ I} are independent.

This property allows to define the transition rate matrix (or infinitesimal generator matrix) Q as:

∀s, s ′ ∈ S, s ̸ = s ′ , Q(s, s ′ ) = P (s, s ′ )E(s) ∀s ∈ S, Q(s, s) = - s ′ ̸ =s P (s, s ′ )E(s) (2.2)
In fact, Property 2.1 is sufficient for the definition of a CTMC.

Proposition 2.1.4 (Infinitesimal definition of a time-homogeneous CTMC)

A stochastic process (S t ) t∈R ≥0 with state space S = {(s

(i) ) i } ⊆ N d is a time-homogenenous
CTMC if and only if there exists a matrix Q = (q ij ) ij called transition rate matrix so that:

∀t ∈ R ≥0 , h > 0, s (i) , s (j) ∈ S, P(S t+h = s (j) |S t = s (i) ) = δ ij + q ij h + o(h).
Remark 2.1.2 (Description of a CTMC.)

A CTMC is fully defined by its initial distribution S 0 and transition rate matrix Q.

E(s) = s ′ ̸ =s Q(s, s ′ )
is the total exit rate of the state s.

The transition rate matrix Q is involved in the evolution of the chain's state over time. Let us define the transient probability vector.

Definition 2.1.13 (Transient state probability)

Considering the countable state space S = {(s (i) ) i }, we define the transient state probability vector of a CTMC as the vector P (t) = (P (t)

i ) i where

P (t) i = P(S t = s (i) ).
The dynamics of the Continuous-Time Markov Chain is described by Kolmogorov's forward equations.

Proposition 2.1.5 (Kolmogorov's forward equations)

Let P be the transient state probability vector (Definition 2.1.13). Then P verifies the Kolmogorov forward equation:

d dt P (t) = P (t) Q.
This equivalence will be proven later in the context of Chemical Reaction Networks (Section 2.3.3).

Remark 2.1.3

For a finite state space (Card(S) < ∞), the solution of the Kolmogorov's forward 2.1. Different perspectives on Markov Chains equations is

P (t) = e Qt .
The analytical form in 2.1.3 is convenient for describing the dynamics of a CTMC system with few states, but the more the number of states increases, the less the computation of the matrix exponential is feasible. This is known as the state-space explosion problem, which justifies the use of statistical and approximate methods for CTMC. The bypass of the state-space explosion problem is the driving force of our work.

CTMC as an oriented graph

In this thesis, we are interested in verification of CTMC, which aims at verifying • Q : S × S → R is the transition rate matrix where ∀s, s ′ ∈S Q(s, s ′ ) = 0.

Remark 2.1.4

In our applications, unless otherwise stated, the initial distribution is combined with a fixed initial state s 0 , which means α = δ s 0 . This case leads to a tuple M = (S, s 0 , Q).

In this case, a CTMC is illustrated by an oriented graph whose nodes are states and edges/transitions represent the transition rate matrix Q. It describes the sleep cycle of a human. The transition rate matrix is defined by

Q(Sleep, Awake) = 1 8 , Q(Awake, Sleep) = 1 16 .
As each row of Q should sum to zero, the two others coefficients equal Q(Sleep, Sleep) = -1 8 and Q(Awake, Awake) = -

Probability measure of CTMCs.

In the following, we describe the set of trajectories of a CMTC and construct a probability measure over this set. It is a requirement for verification of CTMCs: the satisfaction of a formula is based on the computation of this probability measure. We need to ensure that this probability measure is well defined.

Paths/Trajectories of a CTMC.

Definition 2.2.1 (Path/Trajectory of a CTMC)

A path (or trajectory) σ of a CTMC M = (S, α, Q) is a (possibly infinite) sequence

σ = s 0 t 0 -→ . . . t k-1 ---→ s k . . ., with t i ∈ R ≥0 being the sojourn time in state s i ∈ S such that α(s 0 ) > 0 and ∀i, Q(s i , s i+1 ) > 0.
P ath(M) is the set of possible right-continuous paths of a CTMC M.

Several operators are associated with a path σ ∈ P ath(M). For i ∈ N and

t ∈ R ≥0 , σ[i] = s i is the i-th state of σ, δ(σ, i) = t i the sojourn time of σ in the
i-th state, σ@t the state of σ at time t and σ[t⟩ the truncated path of σ at time t (σ[t⟩@0.2 = σ@(t + 0.2)).

For example, if σ = s 0 0.25

--→ s 1 0.5 --→ s 2 0.15 --→ s 3 1 -→ . . . we have σ[1] = s 1 , δ(σ, 2) = 0.15, σ@0.1 = s 1 .
A finite path can be considered as an infinite path with t k = +∞.

Remark 2.2.1 (About the notation of the set of paths)

In the literature, the initial distribution is sometimes not included in the definition of a CTMC. In this context, P ath is used for the set of all paths with any possible initial state and P ath(s) for the paths originating in s.

Probability measure over the set of paths.

Let (S t ) t∈R ≥0 a CTMC. We need a probability measure over the set of paths to quantify how probable a path/trajectory is from a CTMC. In statistical inference, such a probability measure is defined by the pushforward measure of P under S .

Let Φ S be the function:

Φ S : (Ω, F) → (S R ≥0 , F S R ≥0 ) ω → S (ω)
with S R ≥0 the set of right-continuous functions f : R ≥0 → S, and F S R ≥0 the σalgebra generated by S R ≥0 . Then S P = P • Φ -1 S defines a probability measure P S on the measurable set (S R ≥0 , F S T ). We have:

P S = S P = P • Φ -1 S If S is a CTMC, for ω ∈ Ω, S (ω) is a trajectory (like σ). If A is a subset of trajectories, P S (A) = P({ω ∈ Ω, S (ω) ∈ A}).
However, in the representation M = (S, α, Q) of a CTMC, the notion of probability does not appear clearly (except the distribution over the initial state). We only define a transition rate matrix Q, which defines the transitions between states and the timing information of a transition. Intuitively, the elements Q(s, s ′ ) represent the velocity at which the transition s → s ′ occurs. For verification purpose of these systems, the probability measure over the set of trajectories is constructed explicitely by induction. It is based on the assumption that the time before an outgoing transition from a non-absorbing state has an exponential time.

In the following, we define a probability measure P r M , based on the CTMC M = (S, α, Q), over the measurable space (P ath(M), F P ath(M) ), where F P ath(M) is the σ-algebra generated by P ath(M).

Let s i , s j ∈ S so that s i ̸ = s j and E(s i ) > 0. The probability that s i → s j occurs within t ∈ T is given by:

P (s i , s j , t) = P (s i , s j )(1 -e -E(s i )t ) (2.3) where P (s i , s j ) = Q(s i ,s j ) E(s i )
is called a transition probability and (1 -e -E(s i )t ) is the probability that the transition occurs within t.

We define the cylinder set C(s 0 , I 0 , . . . , I k-1 , s k ) with I j =]a j , b j [ as:

C(s 0 , I 0 , . . . , I k-1 , s k ) = {σ ∈ P ath(M), ∀i ∈ {0, . . . , k}, σ[i] = s i , δ(σ, i) ∈ I i } If σ ∈ C(s 0 , I 0 , . . . , I k-1 , s k )
, the sequence of states s 0 → . . . → s k that a trajectory traverses is fixed but the times spent in each state can vary. A way of visualising

C(s 0 , I 0 , . . . , I k-1 , s k ) is s 0 I 0 -→ . . . I k-1 ---→ s k .
The probability measure P r M is defined by induction: P r M (C(s 0 )) =α(s 0 ) = P (s 0 ) P r M (C(s 0 , I 0 , . . . , I k-1 , s k )) =P r M (C(s 0 , I 0 , . . . , I k-2 , s k-1 ))

• P (s k-1 , s k )(e -a k E(s k-1 ) -e -b k E(s k-1 ) ) By denoting T s k ∼ Exp(E(s k )), P r M (C(s 0 , I 0 , . . . , I k-1 , s k )) = P r M (C(s 0 )) • k i=1 P (s i , s i+1 )P(T s i ∈ I i )
The probability measure is thus defined over a subset of F P ath(M) . It is sufficient to define the measure over these cylinder sets: the measure is naturally extended over the whole space (see Appendix B for details).

Proposition 2.2.1 (Consistency of the probability measures) ∀A ∈ F P ath(M) , P r M (A) = P S (A) 2.2. Probability measure of CTMCs.

Proof. It is sufficient to prove the equality of the measures over any cylinder set C(s 0 , I 0 , . . . , I k-1 , s k ). This is a consequence of Caratheodory's extension theorem [START_REF] Klenke | Probability Theory[END_REF], Chapter 1). We refer the readers interested in the probabilistic details of measure constructions to the Appendix B.

Let A = C(s 0 , I 0 , . . . , I k-1 , s k ).
Let us prove this by induction over p ∈ {0, . . . , k}. Let p ≥ 1, A p = C(s 0 , I 0 , . . . , I p-1 , s p ).

Then A k = A.
We set these events:

• E 0 = {S 0 = s 0 } • E 1 = {T s 0 ∈ I 0 , S Ts 0 = s 1 } • j ≥ 2, E j = {a j-1 + j-2 i=0 T s i ≤ T s j-1 + j-2 i=0 T s i ≤ b j-1 + j-2 i=0 T s i , S∑ j-1 i=0 Ts i = s j } = {a j-1 ≤ T s j-1 ≤ b j-1 , S∑ j-1 i=0 Ts i = s j }
E j is the event that after spending a time T s j-1 ∈ I j-1 in s j-1 , we have a jump from s j-1 to s j .

By time-homogeneity of the Markov Chain, P S (A p ) = P(E 0 , . . . , E p ).

Let us reason by finite induction on p ∈ {1, . . . , k}. p = 1:

P(A 1 ) = P(E 1 |E 0 )P(E 0 ) = P(T s 0 ∈ I 0 , S Ts 0 = s 1 |S 0 = s 0 )P(S 0 = s 0 ) = P(T s 0 ∈ I 0 )P(S Ts 0 = s 1 |S 0 = s 0 ) = P (s 0 , s 1 )P(T s 0 ∈ I 0 )
p + 1:

P S (A p+1 ) = P(E 0 , . . . , E p+1 ) = P(E p+1 |E 0 , . . . , E p )P(E 0 , . . . , E p ) = P(E p+1 |E 0 , . . . , E p ) p i=1 P (s i , s i+1 )P(T i ∈ I i )
Then,

P(E p+1 |E 0 , . . . , E p ) = P(T sp ∈ I p , S ∑ p i=0 Ts i = s p+1 |E 0 , . . . , E p ) = P(T sp ∈ I p , S ∑ p i=0 Ts i = s p+1 |S 0 = s 0 , T s 0 ∈ I 0 , . . . , T s p-1 ∈ I p-1 , S∑ p-1 i=0 Ts i = s p )
Since the value of the next state only depends on the previous value, we have:

P(E p+1 |E 0 , . . . , E p ) = P(T sp ∈ I p , S ∑ p i=0 Ts i = s p+1 |S∑ p-1 i=0 Ts i = s p )
Hence, as the CTMC is time-homogeneous, we translate the time byp-1 i=0 T s i , and we obtain by the independence of {T s 0 ∈ I} and {S Ts 0 = s ′ } (Proposition 2.1.3):

P(E p+1 |E 0 , . . . , E p ) = P (s p , s p+1 )P(T sp ∈ I p )
By induction the result is proven.

Remark 2.2.2 (P r M or P S ?)

If these notations appoints the same probability measure, should we drop one of them?

In this thesis, we think both notations are useful. In model checking literature, the representation of Section 2.1.2 and the construction of P r M are frequently presented (e.g. [START_REF] Aziz | Model Checking Continuous Time Markov Chains[END_REF] or [START_REF] Baier | Model-Checking Algorithms for Continuous-Time Markov Chains[END_REF]). For a CTMC M = (S, α, Q), P r M is often denoted P r α (the transition rate matrix is implicit due to the fact we only consider one instance of CTMC in classical model checking).

However, it is convenient to use P S in a statistical context. Also, P r M does not allow probability calculus with random variables, which deprives us of a useful tool.

For example, if we denote A the set of trajectories so that the value at time t 1 is above

x 1 , one can rewrite P r M (A) as P S (A) = P(S t 1 ≥ x 1 ).

Chemical Reaction Networks

In this section, we describe a specific formalism called Chemical Reaction Networks.

They are used to represent a class of biological models that describe the kinetics of chemical phenomena.

In this kind of models, we consider an amount of different species molecules in a specific constant volume: these molecules can collide and provoke transformations.

These phenomena are expressed by reactions:

2.3. Chemical Reaction Networks R i : a i1 X 1 + . . . + a id X d reactants → b i1 X 1 + . . . b id X d products
Due to the complexity of the interactions between the species, the modelling of Chemical Reaction Networks is stochastic. Under particular assumptions, the underlying probabilistic model is a Continuous-Time Markov Chain. This stochastic description is especially useful when the numbers of molecules are not significant.

Most models consider molecules that can collide and react, such as gene regulatory networks [START_REF] Karlebach | Modelling and analysis of gene regulatory networks[END_REF] or cell signaling pathways [START_REF] Eungdamrong | Modeling Cell Signaling Networks[END_REF]. But this framework is also helpful in the study of interactions between humans or animals such as the SIR [START_REF] Kermack | A Contribution to the Mathematical Theory of Epidemics[END_REF] or Lokta-Volterra [START_REF] Lotka | The growth of mixed populations: Two species competing for a common food supply[END_REF]) models.

Definition of a Chemical Reaction Network

In the most detailed description of chemical phenomena, one should describe the position and velocity of each molecule in time and fires a reaction when two molecules collide. Stochastic semantics of Chemical Reaction Network are based on two assumptions:

• The system is well-stirred in a constant volume. At any time t, the system evolution is only represented by a state s ∈ N d where each component s[i] is the number of species. Thus, the system evolution can be described by a stochastic process (S t ) t∈R ≥0 .

• Each possible change of the population occurs through a reaction channel R j .

It is defined by two quantities: ν j the stoichiometric vector and η j the kinetic rate.

The stoichiometric vector ν j represents the variation in number for each species when R j occurs. The kinetic rate η j = η j (s, θ) of a reaction channel R j depends both on the state s and some parameters θ. This quantity is defined so that given a state s, the probability that the reaction j occurs within the infinitesimal interval [t, dt[ is η j (s, θ)dt. In other words,

P(R j occurs within [t, t + dt]| population at time t is s) ≈ η j (s, θ)dt (2.4)
The function s → η j (s, θ)dt is called the propensity function (also called the hazard function).

Definition 2.3.1 (Chemical Reaction Network)

A Chemical Reaction Network is fully described by:

• The set of species S of cardinal d.

• A set R 1 , . . . , R m of reaction channels where each R j is characterised by a pair R j : (ν j , η j ) with ν j = [ν 1j , . . . , ν dj ] the stoichiometric vector, representing the variation in number for each species determined by the occurrence of R j , and η j = η j (s, θ) the kinetic rate of the reaction R j .

• A p-dimensional vector of parameters θ = [θ 1 , . . . , θ p ] affecting the kinetic rate of the reaction channels, with θ in Θ ⊂ R p .

The Equation 2.4 recalls the characterisation of a CTMC in Proposition 2.1.4.

Indeed, the assumptions on the Chemical Reaction Network system evolution leads to a Continuous-Time Markov Chain.

Remark 2.3.1 (CTMC description of a Chemical Reaction Network)

The system evolution of a CRN is described by a CTMC (S t ) t∈R ≥0 :

• A discrete state space S ⊆ N d whose elements are vectors s = [s[1], . . . , s[d]] ∈ S
where s[i] is the population, in terms of number of molecules of the i-th species.

• An initial state s 0 .

• A transition rate matrix Q θ defined by ∀j ∈ {1, . . . , M }, Q θ (s, s + ν j ) = η j (s, θ), and 0 for the states that does not have the form s + ν j .

We deduce the transitition probabilities and total exit rates:

• P θ (s t , s t + ν j ) = η j (st,θ) ηsum(st,θ) • E θ (s) = η sum (s t , θ) = M j=1 η j (s t , θ).
A CTMC whose a state represents the number of individuals of several categories/species is also often called a Markov Population Process.

A common assumption about the kinetic rate is the law of mass action. This law states that the kinetic rate is proportional to the multiplication of the numbers of reactant species.

Remark 2.3.2

Unless stated, the considered Chemical Reaction Networks follows the mass-action law. 20 2.3. Chemical Reaction Networks

Example

Let us illustrate the definitions with an example. Let S = {A, B, C, D} be the set of species, N 4 be the state space. The i-th component i ∈ {1, . . . , 4} of a state is the number of molecules of the species i. We consider the following set of reactions with constant reaction rates θ = (k 1 , k 2 , k 3 ):

R 1 : A + C k 1 -→ D (second order reaction) R 2 : A k 2 -→ B (first order reaction) R 3 : ∅ k 3 -→ C (zero order reaction)
Then the associated stoichiometric vectors are:

ν 1 = [ -1 A is a reactant in R 1 , 0, -1, 1 D is a product in R 1 ] ν 2 = [-1, 1, 0, 0] ν 3 = [0, 0, 1, 0].
The mass-action law leads to the kinetic rates:

η 1 (s, θ) = k 1 • |A| • |C| = k 1 • s[1] • s[3] η 2 (s, θ) = k 2 • |A| = k 2 • s[1] η 3 (s, θ) = k 3 .

Different representations of the system evolution

As we saw in the previous section, the formalism of Chemical Reaction Networks leads to a CTMC (S t ) t∈R ≥0 to describe the system evolution. In the following, we describe the equations that govern (S t ) t∈R ≥0 and how, with additional assumptions, approximated equations can be deduced.

Chemical Master Equation

The Chemical Master Equation of a Chemical Reaction Network, also known as the Kolmogorov equation in probability theory, is the set of equations that describes the probability of each state at any time t.

Let s ∈ S a state. We want to describe the evolution of the probability to be in this state over time.

Let the events E (0) I = "no reaction occurs in the interval I", ∀j ∈ {1, . . . , M }, E (j) I = "the reaction j occurs in the interval time I". Let t > 0 be a time and dt > 0 be an infinitesimal time. Then,

P(S t+dt = s) = P(E (0) [t,t+dt] | S t = s) no reaction occurs P(S t = s)+ M j=1 P(E (j) [t,t+dt] | S t = s -ν j )
one of the reactions occurs

P(S t = s-ν j )
Intuitively, this equation states that the probability of being in state s after an infinitesimal time dt can be decomposed into the sum of two probabilities:

• The probability of being in state s at t and no reaction occurs until t + dt.

• The probability of being in another state so that one of the reactions occurs within dt and brings the system in the state s.

Based on the CTMC framework of a CRN 2.3.1,

P(E (0) [t,t+dt] | S t = s) = (1 - M j=1 η j (s, θ)dt) ∀j ∈ {1, . . . , M }, P(E (j) [t,t+dt] | S t = s -ν j ) = η j (s -ν j , θ)dt
Combining the Equation 2.3.3 and lim dt→0 P(S t+dt =s)-P(S t =s) dt

, we obtain the Chemical Master Equation:

d dt P(S t = s) = M j=1 (η j (s -ν j , θ)P(S t = s -ν j ) -η j (s, θ)P(S t = s)) (2.5)

Random Time Change Representation

In the previous section, we focused on an infinitesimal time to derive an exact equation that describes the evolution in time of our system. Another option is to focus on the total amount of reactions that occur within time t. For j ∈ {1, . . . , M }, N j (S 0 , t) is the random variable that gives the number of occurrences of the j-th reaction until time t. Then,

2.3. Chemical Reaction Networks S t = M j=1
N j (S 0 , t)ν j + S 0 .

( [START_REF] Kurtz | Representations of Markov Processes as Multiparameter Time Changes[END_REF] proves:

S t = M j=1 P( t 0 η j (S v , θ)dv)ν j + S 0 .
(2.6)

where P(a) is a Poisson law with mean a. This is the Random Time Change

Representation of the CTMC.

Tau-leap approximation and Chemical Langevin Equation

Equations 2.5 and 2.6 are in practice difficult to solve. An efficient way to obtain approximate equations that leads to faster simulation algorithms is to discretise the time interval by a step τ . Equation 2.6 becomes:

S t+τ = M j=1 P( t+τ t η j (S v , θ)dv)ν j + S t .
This equation is still exact. Two assumptions (Gillespie, 2007) for the interval time τ can lead to good approximations of our equations:

• For any time t and given the state s t , the propensity function does not vary a lot over

[t, t + τ ], i.e η j (S t ′ , θ) ≈ η(s t , θ) over t ′ ∈ [t, t + τ ]. (Leap condition).
• There is a significant amount of reactions that occur within [t, t + τ ].

These assumptions seem contradictory. Indeed, they assume that τ is small enough so that the propensity function does not vary a lot, but it should be long enough to allow several reactions. In practice, these assumptions are valid if we have enough molecules of each species in the studied system.

The first assumption means that if we know in which state the system is at time t, then the propensity functions are constant over a time interval of length τ . We can then estimate by Riemann sum the integral involved in Equation 2.6 for any N ∈ N:

N •τ 0 η j (S v , θ)dv ≈ N -1 i=0 η j (S iτ , θ)τ
which gives the scheme

S t+τ = M j=1 P(η j (S t , θ)τ )ν j + S t .
(2.7) Given the state s t , the expected number of reactions that will occur is η j (s t , θ)τ .

The second assumption leads to an approximation of a Poisson distribution by a Normal distribution (according to the central limit theorem and the fact that a

Poisson distribution is a sum of unit Poisson distributions).

P(η j (s t , θ)τ ) ≈ N (η j (s t , θ)τ, η j (s t , θ)τ ) = η j (s t , θ)τ + η j (s t , θ)τ N (0, 1)
This gives the scheme:

S t+τ = τ M j=1 η j (s t , θ)ν j + M j=1 ν j η j (s t , θ)τ N (0, 1) + S t (2.8)
which can be seen as a discretisation of the continuous stochastic differential equation:

dS t = M j=1 η j (S t , θ)ν j dt + M j=1 ν j η j (s t , θ)dW (j) t (2.9)
where W (j) t

are independent Brownian motions. The Equation 2.9 is called the Chemical Langevin Equation (Gillespie, 2000).

Reaction Rate Equation: a macroscopic deterministic approximation of a

CRN

Initially, chemical kinetics were described by a set of ODE. The implicit assumption is that the system evolves in a chemical equilibrum. For example, in [START_REF] Kermack | A Contribution to the Mathematical Theory of Epidemics[END_REF] the system evolution of the SIR model is described by three ordinary differential equations (ODE). However, a macroscopic and deterministic equation can be derived from the CME 2.5 that corresponds to the classical ODE 2.3. Chemical Reaction Networks description of chemical systems. We have for any

s (i) ∈ S = {(s (i) ) i∈N } d dt P(S t = s (i) ) = M j=1 (η j (s (i) -ν j , θ)P(S t = s (i) -ν j ) -η j (s (i) , θ)P(S t = s (i) )). (2.10)
As S t is a discrete random variable, E[S t ] = i s (i) P(S t = s (i) ). Then by multiplying the previous equation by s (i) and suming over i ∈ N

d dt E[S t = s (i) ] = M j=1 i∈N η j (s (i) -ν j , θ)s (i) P(S t = s (i) -ν j ) - i∈N η j (s (i) , θ)s (i) P(S t = s (i) ) .
(2.11)

Then, as i∈N η j (s (i) -ν j , θ)s (i) P(S t = s (i) -ν j ) = i∈N η j (s (i) , θ)(s (i) +ν j )P(S t = s (i) ) , it gives:

d dt E[S t ] = M j=1 E[η j (S t , θ)]ν j (2.12)
which gives if we neglect the stochastic variations of our system (t → S t is a deterministic function):

dS t dt = M j=1 η j (S t , θ)ν j .
(2.13) . Equation 2.13 is called the Reaction Rate Equation (RRE). This equation is considered a good approximation of the system's dynamics if the population is very large (Gillespie, 2000).

Remark 2.3.3 (RRE with the Chemical Langevin Equation)

If we assume that the stochastic fluctuations of (S t ) t are negligible, i.e. it is a constant process, then by reconsidering the Chemical Langevin Equation 2.9

dS t = M j=1 η j (S t , θ)ν j deterministic part dt + M j=1 ν j η j (s t , θ)dW (j) t stochastic part . (2.14)
We obtain the Reaction Rate Equation 2.13 by neglecting the stochastic part.

Stochastic simulation of a CRN

We have derived useful exact or approximate equations for the description of the CRN dynamics. Thus, one can derive simulation algorithms of paths from the CRN based on these equations.

Stochastic simulation algorithm

Stochastic Simulation Algorithm (SSA) was first formulated in [START_REF] Gillespie | Exact Stochastic Simulation Of Coupled Chemical-Reactions[END_REF]. It is an exact method that simulates statistically correct paths of a CRN.

In Section 2.1.1, we described the distribution of the sojourn time of a state s t , which is Exponential, and the distribution of transition from one state to another.

These events are independent (Proposition 2.1.3). The probability density function of the couple (S Ts t , T st ) is:

p (S Ts t ,Ts t ) (s t + ν j , τ |S t = s t ) = η j (s t , θ) η sum (s t , θ) Categorical distribution η sum (s t , θ)e -ηsum(st,θ)τ Exponential distribution where η sum (s t , θ) = M j=1 η j (s t , θ) (total exit rate of s t ).
Equation 2.3.4 means that given the state s t at a time t, the time before a reaction occurs and the next reaction that will occur are independent, so we can simulate them step by step, which provides the Stochastic Simulation Algorithm 1.

Chemical Reaction Networks

Algorithm 1 Exact Stochastic Simulation Algorithm Require: A CRN, an end time T, an initial state s 0 .

Ensure: Two collections (s i ) i and (t i ) i that represent the simulated path.

t, t 0 ← 0 i ← 1 while t < T do τ ∼ Exp(η sum (s t , θ)) Take j in {1, . . . , M } with probability ( η j (st,θ) ηsum(st,θ) ) j∈{1,...,M } s i ← s i-1 + ν j i ← i + 1 t, t i ← t + τ end while
The Modified Next Reaction Method is another exact stochastic simulation method developed by [START_REF] Anderson | A modified next reaction method for simulating chemical systems with time dependent propensities and delays[END_REF] based on [START_REF] Gibson | Efficient exact stochastic simulation of chemical systems with many species and many channels[END_REF]. This method needs drawing only one random number per simulated reaction instead of two in the SSA 1. The simulation is based on the representation 2.6. T j (t) = t 0 η j (S v , θ)dv is called the internal time and depends on the propensity function. Algorithm 2 is based on the sequential computation of the internal time, see [START_REF] Anderson | A modified next reaction method for simulating chemical systems with time dependent propensities and delays[END_REF] for more details.

Algorithm 2 Modified Next Reaction Method

Require: A CRN, an end time T, an initial state s 0 .

Ensure: Two collections (s i ) i and (t i ) i that represent the simulated path.

t ← 0 For j ∈ {1, . . . , M }, T j ← 0 i ← 1 while t < T do For j ∈ {1, . . . , M }, P j ∼ Exp(1) For j ∈ {1, . . . , M }, ∆t j ← P j -T j η j (s i-1 ,θ) j * ← arg min {1,...,M } ∆t j s i ← s i-1 + ν j * t ← t + ∆t j * T j ← T j + η j (s i-1 , θ)∆t j * s ∼ Exp(1), P j * ← P j * + s i ← i + 1 end while

Tau-leap approximation

Exact simulation algorithms can be computationally expensive. The greater the end time T , the greater the number of occurred reactions. The idea of the explicit Tau-leaping simulation algorithm is to simulate a bunch of reactions with Poisson distributions at each time step. Indeed, with the tau-leap assumptions, the scheme in Equation 2.7 gives a sequential way to simulate the number of occurrences of each reaction with Poisson distributions during a time τ .

Algorithm 3 Explicit Tau-leaping Simulation Algorithm

Require: A CRN, a number of time steps N , an initial state s 0 .

Ensure: Two collections (s i ) i and (t i ) i that represent the simulated path.

i ← 1 for i = 1:N do ∀j ∈ {1, . . . , M }, k j ∼ P(η j (s (i-1)τ , θ)τ ) s i•τ ← M j=1 k j ν j + s (i-1)τ i ← i + 1 end for
More advanced tau-leaping methods with adaptive τ or implicit schemes can be found [START_REF] Rathinam | Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method[END_REF][START_REF] Cao | Avoiding negative populations in explicit Poisson tau-leaping[END_REF]Cao, Gillespie, and Petzold, 2006).

Example of simulations with the SIR model

To apply the concepts and algorithms we have described earlier in the chapter, we consider the classical SIR compartmental model [START_REF] Kermack | A Contribution to the Mathematical Theory of Epidemics[END_REF].

The SIR model describes the spread of an infectious disease among a population of constant size. The state vector is (X S , X I , X R ) ∈ N 3 . Species S represents the susceptible individuals, I the infected and R the recovered ones. The system's dynamics is described by the Chemical Reaction Network 2.15. Reaction R 1 is the infection process: a susceptible person meets an infected person and gets infected. Reaction R 2 is the recovering process: the infected may become immune to the disease. The parameter vector of the model is θ = (k i , k r ). This chemical system is illustrated by a finite-state CTMC graph with an initial state s 0 = (95, 5, 0) in Figure 2.2.

R 1 : S + I k i -→ 2I R 2 : I kr -→ R (2.15) 2.3. Chemical Reaction Networks (95, 5, 0) (94, 6, 0) (95, 4, 1) (93, 7, 0) (94, 5, 1) (95,1,4) (94,2,4) (95,0,5) k i • 95 • 5 k r • 5
We simulate different paths of the SIR system with θ = (0.12/N pop , 0.05) with SSA and Tau-leaping simulations, and we plot the solution of the reaction rate equation for the infected people.

Figure 2.3 shows the different levels of approximation of the Reaction Rate Equation according to population sizes. On the left, the population is far smaller than on the right. We can see the simulated paths diverges more from the ODE solution than on the right plot. On the right plot, the ODE solution seems a much better approximation.

Figure 2.4 shows the different levels of approximations of the stochastic simulations according to the time step τ . On the right, τ is much higher, and the epidemic's peak is shifted to the right. Possibly a lot of reactions occurs between each time step, so the kinetic rate varies within a time step (violation of one of the tau-leap assumptions). Also, the tau leap paths diverge more from the SSA paths.

Summary

Script 1 (Simulation of the SIR model with SSA, Tau-leaping and ODE.)

This simulation is available in the git repository of this thesis at /code/chap1/sim_ sir.jl.

Summary

In this chapter, we described:

• The fundamentals of Markov Chains in both discrete-time and continuous-time with different points of view.

• The construction of a probability measure of a CTMC over the set of trajectories. It will be used in the case of statistical inference and model checking of CTMCs.

• The Chemical Reaction Network formalism, their stochastic dynamics related to Continuous-Time Markov Chains and classical simulation algorithms of trajectories.

Chemical Reaction Networks are our main subject of study to address both statistical inference and verification tasks which are difficult in most cases. But Markov Chains are also the basis of advanced statistical methods such as MCMC the Bayesian method, which is the next chapter's subject.

Chapter 3

Statistical methods

The general point of parametric statistics is to formalise observations data within a probabilistic parametric model to treat, analyse, infer parameters and make predictions from them.

Statistical inference considers probabilistic models parametrised by a vector θ.

It aims at estimating the parameters based on observations data.

Usually, two points of view are opposed in statistical inference. The first one, the frequentist approach, supposes that a true parameter θ true generates the observations.

Then, statisticians construct an estimator θ(X (1) , . . . , X (n) ) of θ true : it is a function of the random observations X (1) , . . . , X (n) that should converge to θ true and have other convenient statistical properties.

The other one, called Bayesian inference, is the one we are the most focused on in this thesis. In this case, we do not suppose that a parameter θ true generates the data. Instead, we define a prior distribution: it describes our first beliefs about the parameters. Then, this distribution is updated according to observations. The so-called posterior distribution represents our beliefs about the parameters after the integration of observations. We are focused on this type of methods for two main reasons:

• The class of models we are working on is complex, and likelihood computation is, in most cases, intractable. Some groups of methods in the Bayesian approach are likelihood-free. They make likelihood computation dispensable in exchange for less precise inference. Approximate Bayesian Computation is one of them.

They have already proven efficiency in the context of Continuous-Time Markov

Chains [START_REF] Warne | Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art[END_REF][START_REF] Alharbi | Bayesian Inference for Continuous Time Markov Chains[END_REF].

• We use Approximate Bayesian Computation methods to create a statistical verification method for time-bounded reachability problems.

The Bayesian framework

In Bayesian methods, several well-known techniques such as accept-reject, Sequential Monte Carlo, MCMC or Importance Sampling are available to sample parameters from the posterior distribution. However, these methods are not always distinct: they are often connected to create more extensive algorithms to improve performance. Sequential Monte Carlo ABC is an illustrative example. It is thus essential to understand which technique refers to which algorithm to understand the more complex ones. In this Chapter, we describe the different bricks of Bayesian inference to detail a Bayesian inference method we were focused on in our work:

Sequential Monte Carlo ABC.

This Chapter is organised as follows. Section 3.1 details the basic notions for Bayesian inference. Section 3.2 provides multiple Monte Carlo algorithms for sampling from a probability density function π computable up to a constant. Section 3.3 details a family of likelihood-free methods called Approximate Bayesian Computation. Section 3.4 describes kernel density estimation, a non-parametric method to obtain a continuous estimation π of π after we have already managed to sample from the density π.

The Bayesian framework

In mathematical statistics, we consider a statistical model (Definition 3.1.1) of a probability space (Ω, F, P), i.e. a collection of probability distributions.

Definition 3.1.1

Let (Ω, F, P) be a probability space. A statistical model is a collection of probability distributions (P θ ) θ∈Θ indexed by a parameter vector θ.

Usually, statistical models are dominated by a measure that implies a probability density function p θ for any θ ∈ Θ. Thus, the likelihood function is defined based on this density. It measures how probable observations y (1) , . . . , y (N ) are drawn from a probability distribution P θ .

Definition 3.1.2

Let (P θ ) θ∈Θ a parametric statistical model with corresponding densties (p θ ) θ∈Θ . Let y (1) , . . . , y (N ) idependant and identically distributed (i.i.d) observations. The likelihood function is defined as:

l(.; y (1) , . . . , y (N ) ) : Θ → R ≥0 θ → p θ (y (1) , . . . , y (N ) ) = N i=1 p θ (y (i) )
In the frequentist inference approach, the observations are supposed drawn from a distribution P θtrue , i.e. we suppose that one parameter vector θ true generates the data. An estimation of θ true is obtained via the Maximum Likelihood Estimator (Definition 3.1.3). It consists of maximising the log-likelihood. When this problem is too complex, other methods are available such as the Expectation-Maximisation algorithm, which ensures a local minimum.

Definition 3.1.3 (Maximum Likelihood Estimator)

Considering a statistical model (p θ ) θ∈Θ and i.i.d observations y (1) , . . . , y (N ) , a maximum likelihood estimator is a parameter θ ∈ Θ that verifies:

l( θ; y (1) , . . . , y (N ) ) = sup θ∈Θ l(θ; y (1) , . . . , y (N ) ).
In the following, we describe the Bayesian inference framework. This other statistical inference approach is the most focused on in this thesis. As the name mentions, it is based on the Bayes formula, that for two events A and B,

P(A|B) = P(B|A)P(A) P(B) .

.

In the Bayesian approach, uncertainty about the parameters is modelled by a random variable ϑ. P ϑ is a probability distribution called prior distribution. This distribution represents the first beliefs about the parameters before integrating the observations.

In order to formalise the notion of "integration of observations", we need to define the distribution of a random variable conditionally to another, which is embodied by a transition kernel.

Definition 3.1.4 (Transition kernel) Let (X , F X ) and (Y, F Y ) be two measurable spaces. The function Q : X × Y → [0, 1]
is a transition kernel if:

3.2. Monte Carlo methods 1. ∀x ∈ X , Q(x, .) is a probability measure over (Y, F Y ) 2. ∀A ∈ F Y , Q(., A) is measurable.
In the case of two random variables X and Y defined on (X , F X ) and (Y, F Y ),

the transition kernel K that verifies ∀A, B ∈ F X ×F Y , P X,Y (A×B) = X Q(x, B)dP X is called the conditional law of Y given X and it is denoted P Y | X . Remark 3.1.1 (Conditional law of a CTMC)
If we consider a collection of CTMC (S θ ) θ∈Θ , (θ, A) → P S θ (A) is a transition kernel and admits densities.

We consider the random variable Y of observations and a prior distribution P ϑ .

We suppose conditional laws P Y |ϑ and P ϑ|Y admit densities. Bayesian inference is built upon the Bayes theorem:

p(θ|y) = Z -1 p(y|θ)p(θ) ∝ p(y|θ)p(θ)
where Z = θ ′ ∈Θ p(y|θ ′ )p(θ ′ )dθ ′ and ∝ means "proportional to". p(.|y) is called the posterior density of the observation y. This distribution shows how the belief about the parameters is updated after observations of the studied phenomenon. Bayesian estimation deals with the integration problem (instead of the optimisation problem for the frequentist approach) because Z has an integral form. We want to estimate the posterior distribution by estimating the constant Z or even avoiding its computation when it is too complex.

Monte Carlo methods

Monte Carlo methods are a class of algorithmic methods based on simulations of the model. They first appeared after WWII (see [START_REF] Metropolis | The Beginning of the Monte Carlo Method[END_REF] for the emergence of the method). Given a probability P * with density π, the idea is to sample x (1) , . . . , x (n) from π to estimate the expectation of a random variable with density π.

However, simulate samples from π when it is not fully known analytically is difficult. Nicholas Metropolis formulates the first Monte Carlo method based on Markov Chains theory [START_REF] Metropolis | Equation of state calculations for fast computing machines[END_REF][START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Application[END_REF] when π is known up to a constant. Due to the maturity of methods and the fast computing power increase of computers, such related Monte Carlo methods are considered standard techniques in the 2000s [START_REF] Robert | Monte Carlo Statistical Methods[END_REF].

Until the end of the Chapter, the unknown target density to estimate is called π.

Definition 3.2.1 (Density of the empirical distribution)

Let X be a random variable with a density π on X . We consider X (1) , . . . , X (N ) N i.i.d. (independent and identically distributed) observations drawn from π with a realisation x (1) , . . . , x (N ) . The density function of the empirical distribution is defined as:

π(x) = 1 N N i=1 δ x (i) (x) By Glivenko-Cantelli theorem, 1 N N i=1 δ x (i) (x)
a.s.

--→ π(x) which means π is an empirical approximation of the targeted density π. This can lead to an approximation of the expectation for any function h:

E π [h(X)] = X h(x)π(x)dx ≈ X h(x) π(x)dx = 1 N N i=1 h(x i )
The main issue of this estimation procedure is to know how to properly simulate from the density π to get enough samples x (i) in a reasonable time so that the density of the empirical distribution is a good enough approximation of the true probability density function π. Description of the simulation algorithms is the subject of the following sections. The knowledge about π is decreasing as we go along the chapter:

• In Section 3.2.1, π is fully known.

• In Sections 3.2.2, 3.2.3, 3.2.4, 3.2.5, π is known up to a constant, i.e. one can compute a function π(x) ∝ π(x).

• In Section 3.3, the target density π is a possibly intractable posterior distribution. It describes a likelihood-free method that samples approximatively the posterior distribution.

These methods are of much importance in Bayesian inference. Indeed, the classical application of these methods is when π is the posterior distribution p(.|y), and π is p(y|.)p(.) because of the Bayes formula p(θ|y) ∝ p(y|θ)p(θ).

Simulation of a density

The simulation of random variables is based on an available generator of a uniform random variable U (0, 1). Then, each simulation is based on the generalised inverse of a cumulative distribution function (cdf).

Proposition 3.2.1 (Probability integral transform)

Let U = U (0, 1) be a uniform random variable over [0, 1]. Let X be a random variable with cumulative distribution function (cdf) F (F (x) = P(X ≤ x)). Then F -1 (U ) has
the distribution F, where:

F -1 (u) = inf{x/F (x) ≥ u} is the generalised inverse of F .
Any random variable is represented by a transformation of a uniform variable.

If X has a cdf F , it suffices to simulate u ∼ U (0, 1) and take x = F -1 (u). This procedure of simulation requires an explicit computation of F -1 .

Accept-reject algorithm

Sometimes the simulation of the target density π is too complex. The idea of the accept-reject algorithm is to use an auxiliary density q(.) to get samples from π.

Let π(x) ∝ π(x) be an unnormalised density function proportional to the target density. If the ratio π(x) q(x) is bounded, the accept-reject algorithm is guaranteed to generate samples from π, i.e. when:

∃M ≤ 0, ∀x, π(x) < M q(x) (3.1)
Algorithm 4 samples from q(.) and accepts the sample with a probability that depends on π. The inequality 3.1 ensures that the sample is drawn from π.

Algorithm 4 Accept-reject algorithm Require: q proposal density Ensure: (x (i) ) 1≤i≤N N samples drawn from π

for i = 1 : N do repeat x ∼ q(.) u ∼ U (0, 1) until u < π(x) M q(x) x (i) ← x end for
One can prove that the probability of acceptance is

∫ π(x)dx M
, which means M should be as small as possible.

Importance sampling

Another way of simulating the target distribution from an auxiliary density q is the important sampling. In this algorithm, no bound condition is required except that the support of q should include the support of π. The idea is to sample N i.i.d. samples

x (1) , . . . , x (N ) from the importance distribution q, and correct the contribution of each sample x (i) in the empirical distribution 3.2.1 by a weight w(x (i) ) called importance weight. The distribution associated with the Importance Sampling for x (1) , . . . , x (N ) drawn from q is:

π IS (x) = 1 N N i=1 π(x (i) ) q(x (i) ) w(x (i) ) δ x (i) (x)
.

This estimator is unbiased for integration estimation, indeed

E π [h(X)] = X h(x)π(x)dx = X h(x) p(x) q(x) q(x)dx = E q [h π q (X ′ )] ≈ 1 N N i=1 h(x (i) ) π(x (i) ) q(x (i) )
As detailed before, sometimes π is only known up to a constant, i.e. π(x) ∝ π(x).

Monte Carlo methods

To bypass this issue, one can compute the weight with π and normalise it by dividing by the sum of all weights N i=1 w(x (i) ).

Sequential Monte Carlo methods

In complex cases, finding a proposal density q that will aim at the crucial regions of the target density π is difficult. Sequential Monte Carlo methods [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF][START_REF] Moral | Sequential Monte Carlo samplers[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF] were designed to overcome this issue. The idea is to sample from intermediate distributions (π m ) m∈0:M sequentially with proposals (q m ) m∈0:M so that π M = π. These methods are helpful, especially in two cases:

• State-space models. If for one observation x 0:M , M + 1 points are available, it is natural to set π m as the density of the truncated observation p(x 0:m ), m ≤ M .

As x m+1 ∼ p(.|x m ), finding a sequential scheme based on each observation point is more accessible. In this case, the density dimension increases. The target distribution π is the full density p(x 0:M ).

• Densities with decreasing supports. π may have tiny unknown support compared to the space of observations. A strategy is to define a sequence of densities with decreasing support to find more efficiently the small support region of π.

In this case, the densities dimensions may be fixed. This paradigm is used in the Sequential Monte Carlo version of the Approximate Bayesian Computation method (Del Moral, Doucet, and Jasra, 2012a), which we will discuss later. It is related to rare event analysis and subset simulation.

Sequential Importance Sampling

Sequential Importance Sampling is a specific method of Sequential Monte Carlo techniques to overcome issues related to high-dimensional problems related to Importance Sampling. In this part, we suppose each observation x is decomposable, i.e. x = x 0:M .

As mentioned before, we want to sample from intermediate densities (π m ) m∈0:M , i.e. from π 0 to π M = π. At each step m, we sample from a proposal density q m .

A possible choice of proposal densities [START_REF] Moral | Sequential Monte Carlo samplers[END_REF]) at each step is:

q m (x 0:m ) = qm (x m | x 0:m-1 )q m-1 (x 0:m-1 )
The computation of the weight in Importance Sampling for the intermediate distributions can be rewritten as:

w m (x 0:m ) = π m (x 0:m ) q m (x 0:m ) = π m-1 (x 0:m-1 ) q m-1 (x 0:m-1 ) w m-1 (x 0:m-1 ) π m (x 0:m ) π m-1 (x 0:m-1 )q m (x m | x 0:m-1 )
which gives a sequential scheme for the computation. This is the Sequential Importance Sampling method detailed in Algorithm 5. This scheme is particularly interesting in the case of state-space models because the proposal densities can be chosen so that they verify the Markov property, i.e. [START_REF] Moral | Sequential Monte Carlo samplers[END_REF][START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF].

qm (x m |x

0:m-1 ) = q(x m |x m-1 ) (Del
Algorithm 5 Sequential Importance Sampling Algorithm Require: (q m ) m proposal densities, (π m ) m unnormalised intermediate densities Ensure: (x (i) 0:M ) 1≤i≤N N samples drawn from π for i = 1 : N do x (i) 0 ∼ q 0 w (i) 0 ← π0 (x (i) 0 ) q 0 (x (i) 0 ) end for Normalise (w (i) ) i for m = 1 : M do for i = 1 : N do x (i) m ∼ qm (.|x (i) 0:m-1 ) w (i) m ← w (i) m-1 πm(x (i) 0:m ) πm-1 (x (i) 0:m-1 )qm(x (i) m |x (i) 0:m-1 ) end for Normalise (w (i) m ) i end for

Resampling step

As the number of steps M increases, the estimator's variance associated with (S)IS procedures can increase exponentially even in simple cases [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF]. Resampling is a method to overcome this issue.

During the procedure, some samples can have low weights, but our goal is to focus on high probability regions. These low weighted samples will perturb the whole estimation. The more the number of low weights increases, the more the number of valuable samples decreases. This phenomenon is known as particle degeneracy 3.2. Monte Carlo methods [START_REF] Gordon | Novel approach to nonlinear/nongaussian Bayesian state estimation[END_REF][START_REF] Liu | Sequential monte carlo methods for dynamic systems[END_REF]): after some steps, only a few samples will be valuable in our estimation.

One possible solution is to add a step where the low weighted samples are replaced by high weighted samples according to a distribution parametrised by the weights: this is the resampling step. A common way to run this task is multinomial resampling: the new indices of samples (x

(i) m ) 1:N Ensure: Resampling of (x (i) 0:m ) 1:N j 1:N ∼ M(N, (w (i) m ) 1:N ) for i = 1 : N do x (i) 0:m ← x (j i ) 0:m w (i) 0:m ← 1 N end for
One could think this resampling step should be done at every step m to delete the useless samples directly. In pratice, a low weight sample at a step m can have a higher weight at a step m + 1. There is a trade-off between the occurrences of resampling and stability: this is adaptive resampling (Del Moral, Doucet, and Jasra, 2012b).

Effective sample size [START_REF] Kong | Sequential imputations and Bayesian missing data problems[END_REF]) is a statistic to quantify the discrepancy (variability) between the weights. It is defined as:

ESS((w (i) ) i∈1:N ) = 1 N i=1 (w (i) ) 2
.

ESS takes values between 1 (one weight is equal to one and the others 0, worst case) and N (perfectly balanced weights equal 1 N ). Thus, one can use this statistic to construct a condition that will decide when we should resample or not. A threshold

N T is set (typically N T = N
2 ), and if the ESS is lower than the threshold, we resample with Algorithm 7.

Algorithm 7 Adaptive resampling step with Effective Sample Size

Require: (x (i) 0:m ) 1≤i≤N N samples with associated weights (w

(i) m ) 1:N Ensure: Resampling of (x (i) 0:m ) 1≤i≤N N ESS ← ESS((w (i) m ) 1:N ) if N ESS < N T then j 1:N ∼ M(N, (w (i) m ) 1:N ) for i = 1 : N do x (i) 0:m ← x (j i ) 0:m w (i) 0:m ← 1 N end for end if
This resampling step can be then integrated into generic Sequential Monte Carlo procedures. An example with Sequential Importance Sampling is described in Algorithm 8.

Algorithm 8 Sequential Importance Sampling with Resampling Algorithm

Require: (q m ) m proposal densities, (π m ) m unnormalised intermediate densities Ensure: (x (i) 0:M ) 1≤i≤N N samples drawn from π for i = 1 : N do x (i) 0 ∼ q 0 w (i) 0 ← π0 (x (i) 0 ) q 0 (x (i) 0 ) end for Normalise (w (i) ) i Resampling step for (x (i) 0 ) 1:N , (w (i) 0 ) 1:N as in Algorithm 6 or 7 for m = 1 : M do for i = 1 : N do x (i) m ∼ qm (.|x (i) 0:m-1 ) w (i) m ← w (i) m-1 πm(x (i) 0:m ) πm-1 (x (i) 0:m-1 )qm(x (i) m |x (i) 0:m-1 ) end for Normalise (w (i) ) i Resampling step for (x (i) 0:m ) 1:N , (w (i) 0:m ) 1:N as in Algorithm 6 or 7 end for

Markov Chain Monte Carlo

Markov Chain Monte Carlo methods are another type of Monte Carlo methods that sample from a target density π by knowing π ∝ π. In Section 2.1.1 of Chapter 2, we quickly recalled basis of ergodicity theory quickly for discrete state-space Markov

Monte Carlo methods

Chains. The idea of Markov Chain Monte Carlo methods is to define a Markov Chain whose stationary distribution has density π.

However, MCMC methods aim at estimating distributions with continuous state space X . The Markov property of a Markov Chain (X n ) n∈N with a continuous measurable space (X , F X ) is formulated as:

P(X n+1 ∈ A n+1 |X n = x n , X n-1 , . . . , X n ) = P(X n+1 ∈ A| X n = x n ).
The event {X n = x} may have zero probability, so we define a continuous generalisation of n-step transition probabilities (Definition 2.1.3) with the transition kernel (Definition 3.1.4):

P (n) (x, A) = P(X n ∈ A|X 0 = x).
In this context, a stationary distribution Π of a Markov Chain with a transition kernel P is defined as:

Π(A) = X P (x, A)dΠ(x).
Similarly to Theorem 2.1.2 for discrete state-space Markov Chains, ergodic theorems for Markov Chain's convergence to its stationary distribution exist in the continuous state-space case. It is based on Harris Chains, which redefine more generally ergodicity of Markov Chains (Bladt and Nielsen, 2017, Chapter 7.3).

Construction of Markov Chains with a targeted stationary distribution is based on the balance condition.

Proposition 3.2.2 (Balance condition)

Let (X n ) n∈N be a Markov Chain with a continuous state space X and a transition kernel P . We suppose that for all x ∈ X the probability measure on (X , F X ) P (x, .) has a density p(x, y). Let π be a density. If π and p verify the detailed balance condition (or time reversibility conditon):

∀x, y ∈ X , π(x)p(x, y) = π(y)p(y, x) (3.2)
then π is a density of a stationary measure of the Markov Chain.

The Metropolis-Hasting algorithm [START_REF] Metropolis | Equation of state calculations for fast computing machines[END_REF][START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Application[END_REF] is the first formulation of an MCMC algorithm. The idea is to sample a chain from the Markov Chain defined by:

p(x n , x n+1 ) = α(x n , x n+1 )q(x n , x n+1 ), (3.3) where α(x n , x n+1 ) = min(1, π(x n+1 )q(x n+1 , x n ) π(x n )q(x n , x n+1 ) ) = min(1, π(x n+1 )q(x n+1 , x n ) π(x n )q(x n , x n+1 ) )
and for x n ∈ X , q(x n , .) is a proposal density for the next chain's value.

The transition kernel 3.3 verifies the detailed balance condition 3.2, which ensures the Metropolis-Hasting Algorithm's convergence (Algorithm 9) to the stationary distribution π.

Algorithm 9 Metropolis-Hasting Algorithm

Require: A proposal density q(., .), computable π ∝ π Ensure: (x n ) 0≤n≤N drawn from a Markov Chain with limit stationary distribution π Initialise a value x 0

for n = 0 : N -1 do

x * ∼ q(x n , .) u ∼ U (0, 1)

x n+1 ← x * n if u ≤ α(x n , x * ) then x n+1 ← x * end if end for
Another popular MCMC method is Gibbs sampling [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF]. This sampler is useful when the variable of the targeted distribution is decomposable

into several blocks π(x) = π(x[1], . . . , x[k]). Let x n = (x n [1], . . . , x n [k]). The sampling of x n+1 = (x n+1 [1], . . . , x n+1 [k]) is di- vided into k steps, where the i-th step (i ∈ {1, . . . , k}) is: x n+1 [i] ∼ π(.|x n+1 [1], . . . , x n+1 [i-1], x n [i+1], . . . , x n [k])
3.3. Approximate Bayesian Computation: a likelihood-free method Then each sampling step x n [i] is equivalent to a Metropolis-Hasting algorithm step of acceptance probability α = 1 whose proposal distribution q is:

q(x, (x[1], . . . , x[i-1], y, x[i+1], . . . , x[k])) = π(y|x[-i]) = π(x[1], . . . , x[i-1], y, x[i+1], . . . , x[k]) π(x[-i])
where y has the same dimension of x[i] and

x[-i] = x[1], . . . , x[i-1], x[i+1], . . . , x[k].
Thus, the sampling of x n+1 given x n is equivalent to k Metropolis Hasting steps.

One can prove that the target distribution is π [START_REF] Robert | Monte Carlo Statistical Methods[END_REF], which

gives the Algorithm 10.

Algorithm 10 Gibbs Sampling

Require: A target distribution π Ensure: (x n ) 0:N drawn from a Markov Chain with limit stationary distribution π Initialise a value x 0

for n = 0 : N -1 do

x * ← x n for i = 1 : k do y ∼ π(.|x * [-i])
x * [i] ← y end for

x n+1 ← x * end for
More advanced MCMC methods based on Metropolis-Hasting, Gibbs and Sequential Monte Carlo methods such as particle filtering [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] were developed in the literature. We will not detail these methods here, but we refer the reader to [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF], [START_REF] Gilks | Markov Chain Monte Carlo In Practice[END_REF] for a deeper exploration of MCMC methods.

Approximate Bayesian Computation: a likelihoodfree method

So far, the Monte Carlo methods we presented aim at sampling from a distribution π, supposing it is computable up to a constant, i.e. one can compute π ∝ π.

Approximate Bayesian Computation (ABC) is a Bayesian method for the approximate estimation of the posterior distribution when the likelihood is computationally intractable or too expensive. ABC is not the only likelihood-free method; other algorithms exist, such as pseudo-marginal MCMC [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF].

ABC algorithms have gained popularity over the last decade and are applied for parameter inference or model selection in many modelling fields, including systems biology [START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF][START_REF] Ratmann | Using likelihood-free inference to compare evolutionary dynamics of the protein networks of H. pylori and P. falciparum[END_REF][START_REF] Koutroumpas | Bayesian parameter estimation for the Wnt pathway: An infinite mixture models approach[END_REF][START_REF] Lenive | Inferring extrinsic noise from single-cell gene expression data using approximate Bayesian computation[END_REF] or ecology [START_REF] Beaumont | Approximate Bayesian computation in evolution and ecology[END_REF][START_REF] Fasiolo | Approximate methods for dynamic ecological models[END_REF]. The first fundamental work was developed in the field of population genetics [START_REF] Pritchard | Population Growth of Human Y Chromosomes: A Study of Y Chromosome Microsatellites[END_REF], which was motivated by the high dimension and structure of the studied models. ABC methods have proved influential in many applications when classical Bayesian parameter inference methods are challenging to implement, including CTMCs [START_REF] Warne | Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art[END_REF][START_REF] Alharbi | Bayesian Inference for Continuous Time Markov Chains[END_REF].

In this section, we suppose that we observe y exp = (y (1) , . . . , y (N ) ) N i.i.d. observations with y (i) ∈ Y (i) potentially multidimensional, and we consider and prior distribution with density θ → p(θ) and a likelihood p(.|θ). According to Bayes theorem, p(θ|y exp ) ∝ p(y exp |θ)p(θ).

In complex models, p(y exp |θ)p(θ) is not even computable. For example, in the case of CTMCs, we will see that the likelihood for time-discrete observations is intractable in most cases (Section 5.1.2).

ABC Rejection algorithm

The ABC method relies on the same principle as Rejection Algorithm 4 [START_REF] Sisson | Overview of Approximate Bayesian Computation[END_REF], except that it aims at drawing samples from an auxiliary distribution π ABC . The idea consists in drawing parameters from a proposal density q, then simulate the model y ∼ p(.|θ) and accept (θ, y) with a probability proportional to the discrepancy between y and y exp : the closer y is to y exp , the more probable it is to accept the sample.

The quantification of the discrepancy between two samples is based on three functions:

• The summary statistics function η : Y → H. It computes statistics over the simulation result to reduce the dimension of the simulation while maximising the information contained in it.

Approximate Bayesian Computation: a likelihood-free method

• The distance function ρ : H × H → R ≥0 . It is a distance over the summary statistics space H.

• A kernel function K ϵ with a scale parameter ϵ.

Algorithm 11 ABC Rejection Sampling Algorithm Require: q proposal density, K ϵ kernel function with a scale parameter ϵ, M ≥

K ϵ (0) max θ π(θ) q(θ) Ensure: (θ (i) , y (i) ) 1≤i≤n drawn from the ABC posterior π ABC for i = 1 : n do repeat θ ∼ q y ∼ p(.|θ) u ∼ U (0, 1) until u < K (ρ(η(y),η(yexp)))p(θ) M q(θ)
θ (i) , y (i) ← θ, y end for By denoting q(θ, y) = q(θ)p(y|θ) and noticing that:

K ϵ (ρ(η(y), η(y exp )))p(θ) M q(θ) = K ϵ (ρ(η(y), η(y exp )))p(y|θ)p(θ) M q(θ)p(y|θ) ,
Algorithm 11 is equivalent to Algorithm 4 with the proposal density q. Thus, the samples of Algorithm 11 are drawn from the distribution:

π ϵ ABC (θ, y| y exp ) ∝ K ϵ (ρ(η(y), η(y exp )))p(y|θ)p(θ).
It is not clear that this density is an approximation of the posterior distribution.

Let us compute the marginal distribution of θ when η is the identity function. In most of our applications, we use the prior distribution as the proposal density, i.e.

π ϵ ABC (θ|y exp ) ∝ Y π ϵ ABC (θ,
q(θ) = p(θ). Also, the kernel used is the indicator function

K ϵ (u) = 1 2ϵ 1(u ≤ ϵ).
Then, the acceptance probability condition becomes ρ(η(y), η(y exp )). These assumptions lead to a simpler ABC rejection algorithm given in Algorithm 12. In fine, ABC methods aim to draw samples from an auxiliary posterior distribution that approximates the true posterior distribution. Thus, one can use Monte 3.3. Approximate Bayesian Computation: a likelihood-free method Carlo methods developed in 3.2.4 and 3.2.5 with the auxiliary ABC posterior π ϵ ABC as the target density: it is the goal of the following sections.

Markov Chain Monte Carlo ABC

A first improvement of the ABC rejection was proposed by [START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF].

The idea is to create a sampler based on MCMC methods with stationary distribution

π ϵ ABC .
Algorithm 13 Monte Carlo Markov Chain ABC Require: q a proposal density, ϵ, ρ, η Ensure: (θ n ) 0≤n≤N drawn from a Markov Chain with limit stationary distribution

π ϵ ABC .
Find a value θ 0 with Algorithm 12.

for n = 0 : N -1 do θ * ∼ q(.|θ n ) y * ∼ p(.|θ * ) u ∼ U (0, 1) θ n+1 ← θ n if u ≤ p(θ * )q(θn|θ * ) p(θn)q(θ * |θn) and ρ(η(y * ), η(y exp )) ≤ ϵ then θ n+1 ← θ n end if end for

Sequential Monte Carlo ABC

If the tolerance ϵ is tiny, the ABC rejection algorithm can be computationally prohibitive. Thus, an alternative is to create a sampler based on Sequential Monte Carlo In the following, we present the SMC sampler (Algorithm 14) of [START_REF] Beaumont | Adaptive approximate Bayesian computation[END_REF], which is based on [START_REF] Sisson | Sequential Monte Carlo without likelihoods[END_REF]

called Population Monte

Carlo ABC (ABC-PMC). At first, we try to sample from π ϵ 0 ABC with Algorithm 12. ϵ 0 should be big enough to have a moderate initialisation time (otherwise the SMC algorithm is useless). Sometimes, ϵ 0 equals +∞, which is equivalent to sample from the prior. Then, at each step m, a particle θ m is a sample from π ϵm ABC : the parameter is selected. Also, a multinomial resampling procedure is taken into account (Algorithm 6). 50 3.3. Approximate Bayesian Computation: a likelihood-free method Algorithm 14 Population Monte Carlo Algorithm ABC Require: N : number of particles, y exp , (ϵ m ) 0:M , ρ, η, (q m ) 0:M Ensure: (w

(i) M , θ (i) M ) 1≤i≤N weighted samples drawn from π ϵ M ABC Iteration m = 0: find (θ (i) 0 ) 1≤i≤N with Algorithm 12 and ϵ 0 . w (i) 0 ← 1 N , i ∈ {1, . . . , N } for m = 1 : M do for i = 1 : N do repeat Take θ ′ from (θ (j) m-1 ) 1≤jN with probabilities (w (j) m-1 ) 1≤j≤N θ (i) m ∼ qm (.|θ ′ ) y ′ ∼ p(.|θ (i) m ) until ρ(η(y ′ ), η(y exp )) ≤ ϵ m w (i) m ← p ( θ (i) m ) N Σ i ′ =1 w (i ′ ) m-1 qm(θ (i) m |θ (i ′ ) m-1 )
end for

Normalise (w Moral, Doucet, and Jasra, 2012a) also proposes a generic ABC-SMC sampler based on (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], with the possibility of several simulations per particle, an ESS resampling step 7, and an approximation of the weights based on the kernel used in the ABC posterior. However, in the case of a uniform kernel in the acceptance condition, many weights can equal zero after one step (because they do not satisfy the acceptance condition). For a deeper review of ABC samplers, we refer the reader to [START_REF] Fan | ABC Samplers[END_REF][START_REF] Marin | Approximate Bayesian Computational methods[END_REF].

(i) m ) 1≤i≤N end for (Del

Hyperparameters of ABC methods

Summary statistics

When the dimension of y exp is high, to satisfy the acceptance condition is challenging: this issue is known as the curse of dimensionality. Thus, a summary statistics function η is introduced to compute statistics over the observations to reduce the dimension and keep the most complete amount of information contained in the observations. The perfect case is when the statistics are sufficient (Definition 3.3.1). It means that the reduction of the dimension by the summary statistics causes no lack of information.

Unfortunately, the Pitman-Koopman-Darmois theorem states that sufficient statistics whose dimension is bounded while sample size increases only exist for distributions of the exponential family. A practical analytical example with a sufficient statistic is detailed in Appendix C.

Definition 3.3.1 (Sufficient statistic)

Let Y and ϑ be random variables. η is a sufficient statistic if P Y | η(Y )=η(y),ϑ=θ does not depend on θ for any y, θ.

Thus, selecting efficient statistics for the observations is challenging and depends on the type of model knowledge (e.g. for time series [START_REF] Jasra | Approximate Bayesian computation for a class of time series models[END_REF]). It can be seen as a task of information theory: how can we retain maximally the information contained in the samples with a fixed dimension? Thus, several methods based on mutual information, projection techniques, subset selection or auxiliary likelihood were proposed. We refer the reader to [START_REF] Blum | A comparative review of dimension reduction methods in approximate Bayesian computation[END_REF][START_REF] Prangle | Summary Statistics in Approximate Bayesian Computation[END_REF] for a complete review of these methods.

Distance function

The function ρ is a distance function over the summary statistics space H (H = Y when η = id), i.e. it verifies symmetry, separation and triangle inequality. The most common choice is the Euclidean distance:

ρ : H × H → R ≥0 (η (1) , η (2) ) → j (η (1) j -η (2) j ) 2
Recently, considering the simulations as empirical distributions and using distances related to distributions like Kullback-Leibler divergence, Wasserstein or Slice-Wasserstein distances were proposed [START_REF] Bernton | Approximate Bayesian computation with the Wasserstein distance[END_REF][START_REF] Nadjahi | Approximate Bayesian Computation with the Sliced-Wasserstein Distance[END_REF]. In the last two chapters, we will propose new distances to address problems of inference and verification of CTMCs.

Perturbation kernel in ABC-PMC Algorithm

In the SMC sampler of Algorithm 14, a perturbation kernel qm is involved in the local moves of the particles at each step. In most of our computations, we use the multivariate Gaussian kernel:

qm (.|θ ′ ) ∼ N (θ ′ , Σ m )
3.3. Approximate Bayesian Computation: a likelihood-free method where Σ m is the empirical covariance matrix of the parameters (θ (i) m ) i . Other kernels were proposed in [START_REF] Filippi | On optimality of kernels for approximate Bayesian computation using sequential Monte Carlo[END_REF], for example, a multivariate Gaussian kernel with the M nearest neighbours. We tested this method. However, even if it can result in slightly fewer simulations, in practice we observed that the computation of a tree that computes the nearest neighbors at each step increases the algorithm's execution time.

Tolerance level

There are several issues about the tolerance level choice. First, in the SMC version of ABC, a fixed sequence of tolerances (ϵ m ) m is involved, but it is difficult to know in advance which sequence will be computationally efficient. One can set the next ϵ during the execution by considering the α-quantile of the simulation distances from y exp , which gives the adaptive tolerance schedule:

ϵ m+1 = quantile(α, (ρ(η(y (i) ), η(y exp ))) 1≤i≤N ).
In this case, a new hyperparameter α is involved. The difficulty is to choose α so that we do not introduce too many useless steps with a big α (close to 1), but also it should not be too small; otherwise, each step m will be too computationally expensive. Another method developed in (Del Moral, Doucet, and Jasra, 2012a) is to find the new ϵ based on the ESS statistic 7. In their SMC sampler, the weights

w (i)
m depend on ϵ m . Then, the idea is, for a given reduction factor α (typically 0.5 or 0.75),to find the next ϵ by solving the problem:

ESS((w (i) m (ϵ m )) 1≤i≤N ) = αESS((w (i) m-1 (ϵ m-1 )) 1≤i≤N ) over ϵ m ∈]0, ϵ m-1 ].
Thus, we control the particle degeneracy (Section 3.2.4) of our samples with the tolerance evolution.

Also, which final ϵ should we use? We know the smaller the tolerance is, the closer we are to the posterior distribution. A typical practical technique is to fix the running time or the number of simulation and reduce the ϵ as much as we can. However, to our knowledge, no other automatic method exists. Work about asymptotics ABC exists, but it is more related to the concentration of the ABC posterior to a true parameter, which can be seen as a particular case of the posterior distribution when p(.|y exp ) = δ θtrue (.) [START_REF] Frazier | Asymptotic Properties of Approximate Bayesian Computation[END_REF][START_REF] Li | On the Asymptotic Efficiency of Approximate Bayesian Computation Estimators[END_REF][START_REF] Robert | Asymptotics of ABC[END_REF].

Kernel Density Estimation

On the previous sections, we have described algorithms that sample from a target density π, but they do not provide a continuous approximation of π. The goal of kernel density estimation is to get an estimate π with the help of kernel functions.

Indeed, if the empirical distribution 3.2.1 is taken as an approximation, the density π(x) equals zero in most cases, except if x is one of the samples. The key idea of this non-parametric method is that each observation will give mass to the density, i.e. we will construct a density where points near the observations will have positive density.

This mass is given with a kernel function defined below.

Definition 3.4.1 (Kernel function)

A function K : R → R ≥0 is a kernel function if: 1. R K(u)du = 1 2. ∀u ∈ R, K(u) = K(-u)
This section does not intend to present all the existing methods of kernel density estimation. For an exhaustive exploration, we refer the reader to the books (Silverman, 1986) and [START_REF] Chacón | Multivariate Kernel Smoothing and its Applications[END_REF].

Kernel density estimator

We give below the definition of a kernel density estimator [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF].

Definition 3.4.2 (Kernel density estimator)

Let X (1) , . . . , X (N ) be N i.i.d samples from an unknown density π on X . The kernel density estimator π associated with a kernel function K on X is:

∀x ∈ X , π h (x) = 1 N N i=1 K h (x, X (i) )
K h is a rescaled function based on kernel K. The scale factor h is called the bandwidth parameter. h is either a scalar, a vector or a matrix depending on the dimension of X and the choice of K.

In this estimator, each observation gives mass probability on the estimator over the whole set X . The further from the observations, the lower the probability. This definition allows any form for K h . In our work, we have focused on a specific type of kernel density estimator when the density is multidimensional. The product kernel estimator consists in considering the term K h (x, X (i) ) as a product of univariate kernel functions for each dimension of X (i) . It was first defined in [START_REF] Bouezmarni | Nonparametric density estimation for multivariate bounded data[END_REF].

Definition 3.4.3 (Product kernel multivariate density estimator)

Let X (1) , . . . , X (N ) N i.i.d samples from an unknown multidimensional density π on

X ⊂ R d . The s-th component of X (i) is denoted X (i)
s . The product kernel density estimator π associated with kernel functions (K s ) 1≤s≤d on X is:

∀x = (x 1 , . . . , x d ) ∈ X , π(x) = 1 N N i=1 d s=1 K s hs (x s , X (i) s )
where K s is a unidimensional kernel function, K s hs (x s , X

(i) s ) = K s ( xs-X (i) s hs
) and

h = (h 1 , . . . , h s ) ∈ R d ≥0 is the bandwidth vector.
This estimator is flexible and has convenient properties. [START_REF] Bouezmarni | Nonparametric density estimation for multivariate bounded data[END_REF] prove that it is free of boundary bias for specific kernels K s , i.e. the bounds' bias is the same as for interior points of X .

Several choices for the kernel function K are possible and depends on the nature of the estimation problem:

• When X is unbounded, one can preferably use Gaussian kernels 3.4.3,

• When X is bounded with positive probability on bounds, one can preferably use beta or gamma kernels 3.4.3.

Chapter 3. Statistical methods

55

Note that making a supposition on which category our problem belongs to is possible with an histogram visualisation.

Bandwidth selection and Least Squares Cross-Validation

Choosing h is critical to get a good estimator of the density. For the choice of the optimal bandwidth, we consider the minimisation of the mean integrated square error [START_REF] Rosenblatt | Remarks on Some Nonparametric Estimate of a Density Function[END_REF]:

M ISE( π h ) = E π [ R ( π h (x) -π(x)) 2 dx]
This statistic is helpful to prevent our estimator from the bias-variance trade-off.

Analytical computations are often made with order-2 Taylor expansion of the MISE to find the optimal bandwidth that minimises it. For example, if one assumes π is a Gaussian density (π ∼ N (µ, σ)), then the optimal bandwidth is h opt ≈ 1.06σn -1 .

However, in most cases, the analytical expression of the optimal bandwidth is based on the unknown density π, which makes the computation intractable.

Thus, a popular automatic estimation method of the optimal bandwidth is the minimisation of Least Squares Cross-Validation (LSCV). We want to estimate the part of the integrated square error that depends on the bandwidth parameter h.

R( π h ) = X π 2 h -2 X π h π
The first term is the squared L 2 -norm of π h . The second term is classically estimated by Leave-One-Out Cross-Validation [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF].

Definition 3.4.4 (Leave-one out estimator)

The Leave-One-Out estimator of 2 π h π is:

LOO(h) = 2 N N i=1 π (-i) h (X i )
where π

(-i) h (X i ) = 1 N -1 N j=1,j̸ =i K h (X j ).
We can now define the LSCV estimator. The Least Squares Cross-Validation estimator is defined by:

LSCV (h) = X π 2 h -LOO(h)
Thus, the difficulty in the computation of LSCV (h) relies on the first term of the expression. In some cases, such as beta kernels 3.4.3, an analytical expression is not available, and the numerical estimation can be prohibitive depending on the number of observations and the complexity of the kernel. Other methods exist, such as likelihood cross-validation or the test-graph method [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]), but we do not address these techniques in our work.

The optimal bandwidth is chosen by solving:

arg min h LSCV (h)

Kernel functions

As we said previously in the introduction, the choice of the kernel function depends on the nature of the problem, i.e. if X is bounded with positive probability on bounds or not.

Gaussian kernel

For the multivariate Gaussian kernel estimator, we do not use the product of univariate Gaussian kernels, but the multivariate Gaussian kernel. Let X (1) , . . . , X (N ) N i.i.d samples from an unknown density π, X (i) = (X

(i) 1 , . . . , X (i) d ) ∈ X . ∀x ∈ X , π h (x) = 1 N N i=1 1 (2π) d det(H) exp(- 1 2 (x -X (i) ) ⊺ H -1 (x -X (i) ))
where H is a d × d bandwidth matrix.

In this case, the LSCV estimator and its properties are well established [START_REF] Duong | Cross-validation Bandwidth Matrices for Multivariate Kernel Density Estimation[END_REF] and tools such as the R package kernel smoothing (ks) [START_REF] Duong | ks: Kernel Smoothing[END_REF] to perform the estimation.

Beta kernels

In our procedures, we will face density estimation with bounded support. In the following, we focus on the kernel density estimator using beta kernels. They proved their efficiency for estimating probability densities on bounded support when the density has positive probabilities on boundaries.

Thus, we have investigated methods for multivariate KDE with beta kernel and the associated analytical forms of LSCV for beta kernels, which are numerically estimated. This estimation has a non-negligible computation cost.

Univariate beta kernel estimator. In our work, we have to deal with densities that have to be estimated on bounded support, where the probability is possibly not zero on the bounds. To overcome this issue, [START_REF] Chen | Beta kernel estimators for density functions[END_REF]) introduced beta kernel estimators.

Let X (1) , . . . , X (N ) be N i.i.d observations from an unknown density π over

X = [a, b]. We denote the scaled observations X(1) , . . . , X(N) , i.e. X(i) = X (i) -a b-a .
The univariate beta kernel estimator is defined as:

∀x ∈ [a, b], π C1,h (x) = 1 N (b -a) N i=1 K β x-a h(b-a) +1, 1-x-a b-a h +1 ( X(i) )
The kernel function K β α,β (.) is the probability density function of a Beta law with parameters α, β, i.e.

K β α,β (t) = t α-1 (1-t) β-1 B(α,β)
where B is the beta function.

In this case, one can express the ISE:

b a π C1,h (x) 2 dx = b a 1 N (b -a) 2    N i=1 X(i) x-a h(b-a) (1 -X(i) ) 1-x-a b-a h B( x-a h(b-a) + 1, 1-x-a b-a h + 1)    2 dx = 1 N (b -a) 2 N i=1 N j=1 b a X(i) x-a h(b-a) (1 -X(i) ) 1-x-a b-a h X(j) x-a h(b-a) (1 -X(j) ) 1-x-a b-a h B( x-a h(b-a) + 1, 1-x-a b-a h + 1) 2 dx I (i,j) 58

Kernel Density Estimation

With an adequate change of variables,

I (i,j) = (b -a) 1 0 X(i) t h (1 -X(i) ) 1-t h X(j) t h (1 -X(j) ) 1-t h B( t h + 1, 1-t h + 1) 2 dt.
Then,

b a π C1,h (x) 2 dx = 1 N 2 (b -a) N k=1 1 0 X(k) 2t h (1 -X(k) ) 2 1-t h B( t h + 1, 1 -t h + 1) 2 dt + 2 N i=1 N j=i+1 1 0 X(i) t h (1 -X(i) ) 1-t h X(j) t h (1 -X(j) ) 1-t h B( t h + 1, 1-t h + 1) 2 dt
In fine, we want to minimise:

LSCV (h) = 1 N 2 (b -a) N k=1 1 0 X(k) 2t h (1 -X(k) ) 2 1-t h B( t h + 1, 1 -t h + 1) 2 dt + 2 N i=1 N j=i+1 1 0 X(i) t h (1 -X(i) ) 1-t h X(j) t h (1 -X(j) ) 1-t h B( t h + 1, 1-t h + 1) 2 dt - 2 N (b -a) N i=1 1 N -1 N j=1,j̸ =i K β x-a h(b-a) +1, 1-x-a b-a h +1
( X(j) ) [START_REF] Chen | Beta kernel estimators for density functions[END_REF] proves that the optimal bandwidth is of order O(N -2

5 ) as n grows to infinity.

Product kernel estimator and beta kernels. Let X (1) , . . . , X (N ) be N i.i.d samples from an unknown density π on the hyperrectangle

X = [a 1 , b 1 ]×. . .×[a d , b d ] with X (i) = (X (i) 1 , . . . , X (i) d ) ∈ R d .
We denote the scaled observations X(1) , . . . , X(N) with X(i) s = X (i) -as bs-as . Then, one can use the product kernel estimator (Definition 3.4.3) for bounded kernel estimation on X with the beta kernel described in the previous section:

π C1,h (x) = 1 N • V ol(C) N i=1 K β xs-as hs(bs-as) +1, 1-xs-as bs-as hs +1 ( X(i) s ) with V ol(C) = d s=1 (b s -a s ).
Thus, one can express the LCSV estimator as above.

Chapter 3. Statistical methods 59 ( X(j) s ) [START_REF] Bouezmarni | Nonparametric density estimation for multivariate bounded data[END_REF] proves that the optimal bandwidth is of order

LSCV (h) = 1 N 2 V ol(C) ( N k=1 1 0 d s=1 X(k) 2t hs s (1 -X(k) s ) 2 1-t hs B( t hs + 1, 1 -t hs + 1) 2 dt + 2 N i=1 N j=i+1 1 0 d s=1 X(i) t hs s (1 -X(i) s ) 1-t hs X(j) t hs s (1 -X(j) s ) 1-t hs B( t hs + 1, 1-t hs + 1) 2 dt) - 2 N • V ol(C) N i=1 1 N -1 N j=1,j̸ =i K β
O(N -2 d+4 ).
Remark 3.4.1 (Gamma kernels)

For estimation problems when X is partially bounded (e.g. X = [0, +∞[ in the univariate case), one can use gamma kernels

(K γ k,θ (t) = t k-1 exp( -x θ ) Γ(k)θ k
where Γ is the gamma function) as proposed in [START_REF] Chen | Probability Density Function Estimation Using Gamma Kernels[END_REF][START_REF] Bouezmarni | Nonparametric density estimation for multivariate bounded data[END_REF].

The methodology is the same as developed above.

Estimation of the bandwidth in practice. The estimation of the bandwidth parameter induces an optimisation problem:

arg min h LSCV (h)
The bigger n, the more computationally expensive LSCV (h) is. Thus, we consider the following gradient-free procedures for product kernel estimators:

• For the univariate beta kernel estimator, we consider Brent's method [START_REF] Brent | An Algorithm with Guaranteed Convergence for Finding a Zero of a Function[END_REF]).

• For the multivariate beta kernel, we consider the simulated annealing algorithm [START_REF] Goffe | Global Optimization of Statistical Functions with Simulated Annealing[END_REF].

• If these methods are too computationally expensive, one can consider a grid search around the convergence rate of the bandwidth O(N -2 d+4 ). This procedure is easily distributable.

Script 2 (Bounded Kernel Density Estimation in Julia language)

These procedures were implemented in Julia, accesible at the git repository https: //gitlab-research.centralesupelec.fr/2017bentrioum/boundedkde.jl with a 3.5. Summary focus on computation efficiency given that no implementation of product kernel estimators was available in Julia.

Remark 3.4.2 (π is supposed unknown)

In the whole section, we have described non-parametric methods to estimate π when the density is unknown. Suppose one has an intuition about the family of densities to which the unknown density π belongs. In that case, one can consider a parametrised family of probability density functions and fit the observations to reduce the complexity of the problem.

Summary

In this chapter, we described:

• The Bayesian framework and most common Monte Carlo methods that sample from an unknown density π when we can compute this density up to a constant.

• A family of likelihood-free methods called Approximate Bayesian Computation, which approximately samples from a posterior density, even if the likelihood is computationally intractable.

• Kernel density estimation methods for the estimation of an unknown density based on its samples, focusing on bounded support distributions.

Chapter 4

Verification of Continuous-Time

Markov Chains

The world contains many complex systems such as communication networks, industries or transport that need to be safe, i.e. verify robust specifications. These systems may be modelled by deterministic or probabilistic models. Model checking [START_REF] Baier | Principles of model checking[END_REF] groups several techniques for formal verification of models. This statement raises two questions:

• How a specification is formally described?

• Which class of models can be considered? How the formulation of model checking problems evolves along with the class of models?

In this Chapter, we focus on the verification of Continuous-Time Markov Chains. This is motivated by the growing interest in model checking in Systems Biology [START_REF] Kwiatkowska | Using Probabilistic Model Checking in Systems Biology[END_REF], which our subject study, Chemical Reaction Network models, belongs to. In particular, we focus on the statistical methods related to these verification tasks. Indeed, in the last two decades, statistical methods have emerged to balance the state-space explosion from which numerical methods of

Temporal logic

To verify the specifications of our systems, we need a logical formalism to describe them. In this section, we describe the syntax and the satisfaction relation of two temporal logics: Metric-Interval Temporal Logic (MITL) [START_REF] Alur | The benefits of relaxing punctuality[END_REF], which is based on the Linear Temporal Logic (LTL) and does not involve probabilistic operator; and CSL [START_REF] Aziz | Model Checking Continuous Time Markov Chains[END_REF][START_REF] Baier | Model-Checking Algorithms for Continuous-Time Markov Chains[END_REF], which is based on the Computation Tree Logic (CTL).

In the following, we consider a CTMC M = (S, α, Q) and its trajectory set P ath(M).

MITL Definition (Metric Interval Temporal Logic)

An MITL formula φ is a term of the following grammar:

φ ::= true | µ | ¬φ | φ 1 ∧ φ 2 | φ 1 U I φ 2
where true stands for the true formula, µ ∈ AP is an atomic proposition, ¬ and ∧ are the negation and conjunction operators of propositional logic, and U I is the until temporal operator with I ⊆ R ≥0 a non-singular interval. The satisfaction relation |= of an MITL temporal formula is defined for a path σ ∈ P ath(M) as follows:

σ |= true σ |= µ ⇔ σ[0] |= µ σ |= ¬φ ⇔ σ ̸ |= φ σ |= φ 1 ∧ φ 2 ⇔ σ |= φ 1 and σ |= φ 2 σ |= φ 1 U I φ 2 ⇔ ∃t ∈ I, σ[t⟩ |= φ 2 and ∀t ′ < t, σ[t ′ ⟩ |= φ 1
For s a state of the CTMC M, s |= φ reads: state s satisfies φ. Although MITL semantics is defined for paths of a CTMC the truth of an MITL formula is naturally extended to states of a CTMC by considering the set of paths originating from a given path (see the second row of Definition 4.1.2).

For example, as

σ[t⟩ |= µ ⇔ σ@t |= µ, a time-bounded until formula φ 1 U [t 1 ,t 2 ] φ 2
is satisfied at time t for a path σ if and only if there exists

t ′ ∈ [t 1 , t 2 ] such that σ@(t + t ′ ) |= φ 2 and ∀t ′′ < t ′ , σ@(t + t ′′ ) |= φ 1 , which is expressed as σ@t |= φ 1 U [t 1 ,t 2 ] φ 2 . σ |= φ means σ@0 |= φ 1 U [t 1 ,t 2 ] φ 2 as stated in Definition 4.1.2.

CSL

Continuous Stochastic Logic (CSL) [START_REF] Aziz | Model Checking Continuous Time Markov Chains[END_REF][START_REF] Baier | Model-Checking Algorithms for Continuous-Time Markov Chains[END_REF] is a logic derived from CTL and adapted to Continuous-Time Markov Chains.

Definition 4.1.3 (Continuous Signal Logic)

A CSL formula ϕ is described by the following grammar:

ϕ ::= true | µ | ¬ϕ | ϕ ∧ ϕ | ϕ ∨ ϕ | P ▷◁p φ | S ▷◁p ϕ φ ::= ϕU I ϕ | X I ϕ
where ▷◁∈ <, ≤, ≥, > and I ⊆ R ≥0 is a non-singular interval.

ϕ is a state formula, whereas φ is a path formula. For path formulae, CSL has an Until operator U, and a next operator X that asserts if the next state of the path within the interval satisfies the property.

CSL includes two operators: the probability operator P and the steady-state operator S (also called long run operator). P ▷◁p φ means the probability that a path verifies φ is ▷◁ p. S ▷◁p ϕ means the steady-state probability that ϕ holds is ▷◁ p.

Model checking of Continuous-Time Markov Chains

The two types of formula lead to the definition of two satisfaction relations, one for state formulae and another for path formulae. The satisfaction relation |= over state formulae is defined as: The satisfaction relation |= over path formulae is defined as:

s |= true s |= ¬ϕ ⇔ s ̸ |= ϕ s |= ϕ 1 ∧ ϕ 2 ⇔ s |= ϕ 1 ∧ s |= ϕ 2 s |= P ▷◁p φ ⇔ P r M 0 ({σ ∈ P ath(M 0 ), σ |= φ}) ▷◁ p with M 0 = (S,
σ |= X I ϕ ⇔ σ[1] |= ϕ ∧ δ(σ, 0) ∈ I σ |= ϕ 1 U I ϕ 2 ⇔ ∃t ∈ I, σ@t |= ϕ 2 ∧ ∀t ′ < t, σ@t ′ |= ϕ 1

Eventually and global operators

The Until operator is less convenient if one wants to express properties like "the variable X reaches x 1 between I" or "the variable X is contained in [x 1 , x 2 ] during I". Then, we consider two derivations of the Until operator: the eventually F I φ = true U I φ, which stands for "at some point within I φ is satisfied", and the globally G I φ = ¬F I ¬φ which stands for "φ is always satisfied within I".

Model checking of Continuous-Time Markov Chains

Now that we have defined how to describe specifications for CTMCs with temporal logic, we detail how we verify these specifications. Considering a specification ϕ, we want to know if a CTMC M fulfils ϕ. It relies on the computation of the probability that a path formula φ is satisfied. In the model checking literature, the threshold problem is sometimes referred to as the qualitative approach; whereas the estimation problem is sometimes referred to as the quantitative approach of model checking.

About numerical methods

Numerical methods are exact verification methods that exhaustively compute the probabilities contained in a CSL formula ϕ over the whole CTMC state space. It implies the evaluation of sub-formulae over the states of the CTMC to label the states with the entire formula ϕ. For example, if ϕ = ¬ϕ 2 , the verification is done by labelling with ϕ the states that are not labelled with ϕ 2 .

For example, in the case of a CSL formula ϕ = P ▷◁p φ with φ a time-bounded Until formula, the verification of ϕ implies transient analysis, i.e. the computation of the transient probabilities described by Kolmogorov's Forward Equations 2.1.5.

Model checking of Continuous-Time Markov Chains

As we have seen before, these methods can suffer from the state-space explosion, adding that in practice, the required memory to run these algorithms is high and the parallelism difficult. Several classical improvements exist for faster model checking as the aggregation of states (reduction of the state space) or uniformisation of CTMC [START_REF] Baier | Model-Checking Algorithms for Continuous-Time Markov Chains[END_REF], which we will not detail here.

Numerical methods for model checking of CTMC are available in the Prism model checker software [START_REF] Kwiatkowska | {PRISM} 4.0: Verification of Probabilistic Real-time Systems[END_REF].

Statistical model checking

Estimation problem -Confidence bounds

We consider the estimation problem in Definition 4.2.2. We want to estimate p = P r(φ; M). Simulating σ ∼ M and verifying if σ |= φ can be seen as a draw from a Bernoulli random variable X ∼ B(p). Indeed, there is a probability p that σ |= φ.

Statistical model checking aims at estimating p with statistical confidence: we want to know how many simulation runs we need to have a certain level of statistical confidence. This task is done by the computation of a confidence interval.

Definition 4.2.4

Let X be a random variable and p a parameter to estimate. A confidence interval of

level (1 -α) is a random interval [b(X), u(X)] with P(b(X) ≤ p ≤ u(X)) ≥ 1 -α. For our Bernoulli variable X ∼ B(p), E[X] = p. A natural estimator is X N = X (1) +...+X (N ) N
which converges to p.

The Hoeffding inequality gives for any c ∈ R ≥0

P(|X

N -X| ≥ c) ≤ 2 exp(-2c 2 N ) P(|X N -X| ≤ c) ≥ 1 -2 exp(-2c 2 N )
By setting c so that 2 exp(-2c 2 N ) = α, we obtain the following confidence interval: of the desired confidence interval. Dividing r by two leads to run four times more simulations. In [START_REF] Jegourel | Sequential schemes for frequentist estimation of properties in statistical model checking[END_REF], they use Massart bounds to sharpen the bounds of the confidence interval.

P X N - -log(α/2) 2N ≤ p ≤ X N + -log(α/2) 2N ≥ 1 -α.
When p is tiny (less than 10 -4 ), Statistical Model Checking is combined with a specific field of statistics called rare event analysis. It deals with methods such as subset simulation and important sampling mentioned in Chapter 3 (e.g. see [START_REF] Barbot | Acceleration for Statistical Model Checking[END_REF]). For a complete review of Statistical Model Checking methods, we refer the reader to [START_REF] Legay | Statistical model checking: An overview[END_REF][START_REF] Legay | Statistical Model Checking[END_REF].

Threshold problem -Hypothesis testing

We consider a CSL formula ϕ = P ≤p φ and the threshold problem of Definition 4.2.3.

If numerical methods are too computationally expensive, a method to statistically verify this specification is hypothesis testing.

The principle of hypothesis testing is to formulate two disjoint hypotheses and test their validity based on observation data X (1) , . . . , X (N ) . The first hypothesis H 0 , called the null hypothesis, is a default position and H 1 the alternative one. Then, we compute a statistic T (X (1) , . . . , X (N ) ), for which we can derive the exact or asymptotic distribution under the null hypothesis. If the realisation of this statistic with our data is too unlikely compared to the distribution under the null hypothesis, the null hypothesis is rejected.

For statistical model checking, the most used statistic is Wald's Sequential Probability Ratio Test. We define the two hypotheses as:

• H 0 : p ≥ p 0 • H 1 : p ≤ p 1
Let x (1) , . . . , x (N ) be a realisation of X (1) , . . . , X (N ) ∼ B(p). Then, the test statistic is the quantity

t = p 1,N p 0,N = N i=1 P(X (i) = x (i) |p = p 1 ) P(X (i) = x (i) |p = p 0 ) = p d N 1 (1 -p 1 ) N -d N p d N 0 (1 -p 0 ) B-d N
where d N = N i=1 x (i) . t is computable sequentially, which is convienent to obtain a decision without knowing the sufficient number of simulations in advance.

Model checking of parametric Continuous-Time Markov Chains

Let α be the probability of the type I error (rejection of the true null hypothesis) and β be the probability of the type II error (not rejecting while H 0 is false). Then, the decision process is:

1. If t ≤ β 1-α , accept H 0 . 2. Else if t ≥ 1-β α , accept H 1 .
3. Else, no conclusion.

Model checking of parametric Continuous-Time Markov

Chains

The previous section's methods deal with an instance of a CTMC. In some cases, Let Θ be a non-empty set. A parametric CTMC, abbreviated pCTMC, (M θ ) θ∈Θ (or

(S θ ) θ∈Θ
) is a collection of CTMCs indexed by a parameter vector θ ∈ Θ.

The general goal of these methods is to analyse how the satisfaction probability of temporal properties evolve along with the parameters.

Estimation problem -Satisfaction function regression

We extend the verification task described in Definition 4.2.2. We no longer want to estimate a value P r(φ; M) for a specific value θ but rather a set of values (P r(φ; M θ )) θ∈Θ .

Definition 4.3.2 (Satisfaction probability function)

Let (M θ ) θ∈Θ be a parametric CTMC and φ an MITL formula. The function

f φ :Θ → [0, 1] θ → P r(φ; M θ )
is called the satisfaction probability function.

The satisfaction probability function (Definition 4.3.2) expresses how the satisfaction probability of φ varies with the parameters. Then, the estimation problem for parametric model checking is to estimate the satisfaction probability function.

If Θ is finite, then the task can be tackled by iterating numerical or statistical methods of classical model checking for each parameter. If Θ is infinite, numerical methods are impossible (it requires an infinite computation time). The problem becomes a task of function regression. An exact method means the analytical computation of transient probabilities is possible for every parameter θ ∈ Θ. In most cases, it is impossible (we give an example that is simple enough to have an analytical solution in Section 6.5.1). Such difficulties encourage statistical methods.

Remark 4.3.1 (Automaton-ABC method)

Chapter 6 details a new method for the estimation of the satisfaction probability function based on ABC and Linear Hybrid Automata.

Statistical methods

In [START_REF] Ceška | Precise Parameter Synthesis for Stochastic Biochemical Systems[END_REF], they formulate the max synthesis, aiming to estimate the regions of Θ where the satisfaction probability function 4.3.2 is maximal, with confidence bounds. A significant work about the estimation of the satisfaction probability function for parametric CTMC 4.3.1 has been done in [START_REF] Bortolussi | Smoothed model checking for uncertain Continuous-Time Markov Chains[END_REF]. The method is called Smoothed Model Checking. Smoothed MC algorithm estimates the satisfaction probability function in a Bayesian framework with the help of Gaussian Processes (GP).

Let us define the method's framework based on Bayesian statistics. We want to estimate the function f φ based on random observations denoted O (1) , . . . , O (n) whose realisations are o (1) , . . . , o (n) (the meaning of these observations is detailed below).

Smoothed model checking is a Bayesian learning regression method. We have to define 1. The prior f , 2. The observation model.

Smoothed MC estimates f φ by estimating the posterior distribution:

p(f (θ * )| O (1) = o (1) , . . . , O (n) = o (n) ), θ * ∈ Θ.

Model checking of parametric Continuous-Time Markov Chains

In function regression, defining a prior amounts to define a stochastic process (also called random function) (f (θ)) θ∈Θ . It is analogous to the prior definition in statistical parameter inference when one defines a random variable over the parameters space. In the Smoothed MC method, the prior is set as a Gaussian process defined below.

Definition 4.3.3 (Gaussian process)

A Gaussian process is a stochastic process (Y x ) x∈X whose finite dimensional distributions are Gaussian, i.e.

∀n ∈ N, ∀x 1 , . . . , x n ∈ X , (Y x 1 , . . . , Y xn ) ∼ N (µ, Σ)
It is uniquely defined by its mean function µ(x) = E[Y x ], and its covariance function

k(x, x ′ ) = Cov(Y x , Y x ′ ).
We can write Y ∼ GP (m(.), k(., .)).

A Gaussian process is entirely described by its mean function and covariance function. The mean function is typically constant (either zero or the mean of the dataset). The chosen covariance function in Smoothed MC is the Gaussian Radial

Basis function k(x, x ′ ) = σ 2 e -∥x-x ′ ∥ 2 /λ 2 .
Second, we consider n observations o (1) , . . . , o (n) related to parameters θ (1) , . . . , θ (n) .

Each o (i) is the number of trajectories that verifies φ among the simulation of M trajectories with the CTMC M θ (i) (as it can be done in statistical model checking, Section 4.2.2). Thus, o (i) is the result of M independent drawing from a Bernoulli law B(f φ (θ (i) )), it is a Binomial law Bin(M, f φ (θ (i) )).

O (i) | f (θ (i) ) ∼ Bin(M, f (θ (i) ))
Thanks to Bayes rule,

p(f (θ * ), f (θ (1) ), . . . , f (θ (n) )| o (1) , . . . , o (n) ) ∝ p(f (θ * ), f (θ (1) ), . . . , f (θ (n) ))
Gaussian process

• n i=1 p(o (i) | f (θ (i) )) Bin(M,f (θ (i) ))
Then, the prediction is performed by the marginalisation of the density written below.

p(f (θ * )| o (1) , . . . , o (n) ) = Θ n p(f (θ * ), f (θ (1) ), . . . , f (θ (n) )| o (1) , . . . , o (n) )df (θ (1) ) . . . df (θ (n) ) (4.1)
which is analytically intractable. Smoothed MC uses the Expectation-Propagation algorithm [START_REF] Minka | Expectation Propagation for approximate Bayesian inference[END_REF] to compute an approximate posterior and control the error.

To

Y x 1:N Y x * ∼ N        0,        (k(Y x i , Y x j )) 1:N k(Y x 1 , Y x * ) . . . k(Y x N , Y x * ) k(Y x * , Y x 1 ) . . . k(Y x * , Y x N ) k(Y x * , Y x * )              
Thus, the conditional law Y x * |Y x 1:N is easily computable with Gaussian identities, and the expectation has a closed form.

Parameter synthesis -threshold problem

Parameter synthesis aims at partitioning the parameter space Θ with a CSL formula ϕ. We want to estimate two subsets of parameters Θ ϕ and Θ ¬ϕ , such that Θ ϕ contains the parameters for which M θ satisfies the specifications, and conversely for Θ ¬ϕ .

•

Θ ϕ = {θ ∈ Θ, M θ |= ϕ} • Θ ¬ϕ = {θ ∈ Θ, M θ ̸ |= ϕ}
which can be seen as a classification task. In practice, it is difficult to solve this problem exactly. In the case of an infinite parameter space Θ, it would require that transient analysis can be solved analytically. Thus, the parameter synthesis (also called threshold synthesis) is formulated as follows.

Let ϕ be a CSL formula let ϵ > 0 be a tolerance error. We want to find a partition

(T, F, U ) of Θ s.t. 4.4. Hybrid Automata Stochastic Logic • T ⊆ Θ ϕ • F ⊆ Θ ¬ϕ • U ⊂ Θ with V olume(U )/V olume(Θ) ≤ ϵ
Several numerical methods have been published to solve this problem. [START_REF] Han | Approximate parameter synthesis for probabilistic time-bounded reachability[END_REF]) address an approximate solution by discretising Θ for time-bounded reachability probabilities. In [START_REF] Ceška | Precise Parameter Synthesis for Stochastic Biochemical Systems[END_REF], they solve the problem by computing the bounds of transient probabilities of [START_REF] Brim | Exploring Parameter Space of Stochastic Biochemical Systems Using Quantitative Model Checking[END_REF]).

Statistical formulation

If numerical methods are too computationally expensive, we can formulate a statistical version of threshold synthesis.

Let ϕ = P ≥p φ with φ a path formula. Considering an estimator f φ (a random function) of the true satisfaction probability function f φ , a tolerance ϵ, a confidence threshold α, we want to find a partition ( T , F , Û ) of Θ s.t.

• T = {θ ∈ Θ, P( f φ (θ) ≥ p) ≥ 1 -α} • F = {θ ∈ Θ, P( f φ (θ) < p) ≥ 1 -α} • Û ⊂ Θ with V olume(U )/V olume(Θ) ≤ ϵ
To our knowledge, a statistical approach of parameter synthesis has only been tackled by [START_REF] Bortolussi | Bayesian statistical parameter synthesis for linear temporal properties of stochastic models[END_REF], which uses the Smoothed MC estimator of the function to compute statistical bounds. In [START_REF] Molyneux | ABC ( SMC ) 2 : Simultaneous Inference and Model Checking of Chemical[END_REF], they use Support-Vector Machines to approximate the partition of Θ according to the CSL formula. Indeed, parameter synthesis is a classification task. Thus any supervised learning method is worth considering.

Hybrid Automata Stochastic Logic

Hybrid Automata Stochastic Logic (HASL) is a logic [START_REF] Ballarini | HASL: An expressive language for statistical verification of stochastic models[END_REF] designed for Discrete-Event Stochastic Processes (DESP), including CTMCs. It is a quantitative logic that allows operators to compute real-valued indicators like expectation or maximum. This logic is more expressive than CSL for non-nested formulae (Ballarini et al., 2011, Section 3.3). This quantitative logic inherently adopts the statistical point of view of verification, and is based on model simulations, which makes it wellsuited for statistical model checking of continuous-time models with discrete events.

An HASL formula is described by a couple (Z, A), where A is a linear hybrid automaton A, and Z is an HASL expression.

Stochastic Petri Net

The implementation of HASL (see Section 4.4.4 further) is based on Generalised Stochastic Petri Nets (GSPN) formalism. Chemical Reaction Network models can be described by Stochastic Petri Nets (SPN), which is a special case of GSPNs.

Definition 4.4.1 (Stochastic Petri Net)

A Stochastic Petri Net is a tuple SP N = (P, T, F, M 0 , Λ) where

• P = {p 1 , . . . , p d } is a set of places.
• T = {t 1 , . . . , t M } is a set of transitions.

• W : (P × T ) ∪ (T × P ) → N is a function that assigns a multiplicity to each possible arc.

• m 0 ∈ N d is the initial marking, which represents the initial distribution of tokens over the places.

• Λ = {λ 1 (.), . . . , λ M (.)} is a vector of firing rates that may depend on the marking of the SPN.

An SPN is a bipartite graph. A place p can only be connected to a transition t and vice versa, which is represented by W (p, t) > 0 (or W (t, p) > 0). A marking m is a vector of dimension d. It assigns an integer to each place. Each integer represents the number of tokens contained in a place. A firing rate λ(m) is a transition rate of some transition t i , i.e. the firing time of t i is distributed as Exp(λ(m)).

An SPN easily describes a Chemical Reaction Network (Section 2.3, Chapter 2).

Indeed, the transitions are the chemical reactions, i.e. T = {R 1 , . . . , R M }. Each place represents a species, the number of tokens contained in a place is the population of species, and a marking m is an element s ∈ S.

For a marking m, the stochastic dynamics of an SPN are composed of two steps:

Hybrid Automata Stochastic Logic

• Enabling of transitions. A transition t i is enabled if all its input places p (which means W (p, t i ) > 0) have a number of tokens greater than W (p, t i ).

• Firing of a transition. A enabled transition is fired if it has the minimum sojourn time within the enabled transitions. The sojourn times are exponentially distributed with mean 1 λ i (m) . This is analogous to the SSA (Algorithm 1) over the enabled transitions.

Let us illustrate this with the SIR Chemical Reaction Network whose equations are recalled below: This SPN has two transitions R 1 and R 2 . The place S has 95 tokens, I has 5 tokens, whereas no tokens are in R. S and I have at least one token, which enables R 1 . Also, I has at least one token, which enables R 2 . If R 1 occurs, one token is moved from the place S (W (S, R 1 ) = 1), whereas the place I receives a token.

R 1 : S + I k i -→ 2I R 2 : I kr -→ R
(W (R 1 , I) -W (I, R 1 ) = 1). The marking is updated to m 1 = (94, 6, 0).

Remark 4.4.1 (Generalised Stochastic Petri Nets)

GSPNs are SPNs that also include immediate transitions. The analogous stochastic processes are Markov Renewal Processes. We refer the interested readers to [START_REF] Bause | Stochastic Petri Nets -An Introduction to the Theory[END_REF] for a deeper exploration of Petri Net models and their timed extensions.

Linear Hybrid Automata

Definition

Linear Hybrid Automata (LHA) are a class of automata that extends Deterministic Timed Automata (DTA). They are equipped with clocks and can deal with continuous variables as well as discrete changes.

Definition 4.4.2 (Linear Hybrid Automaton)

A Linear Hybrid Automaton (LHA) is a tuple A = (E, L, Λ, Init, Final, X, Flow, →)

where

• E is a finite alphabet of events.

• L = {l 0 , . . . , l q } is a finite set of locations.

• Λ : L → P rop is a location labelling function.

• Init ⊂ L is the set of initial locations.

• Final ⊂ L is the set of final locations.

• X = (x 1 , . . . , x n ) is a n-tuple of data variables.

• Flow : L → Ind n is a function which associates each location with n indicator functions (one per data variable). Flow i is the i-th projection.

• The transition relation (the set of edges) →, which is a subset of

L × (2 E × Const) ⊎ ({♯} × lConst) × U p × L where
-♯ means the transition is asynchronous, i.e. the transition can be fired even if it is not synchronised with an event of the stochastic model, -Const is a set of constraints, -lConst is a set of left-closed linear constraints.

A transition (l, (E, γ), U, l ′ ) ∈→ is written l

E,γ,U ----→ l ′ .
This definition is coupled with additional properties.

• Initial determinism: ∀l, l ′ ∈ Init with l ̸ = l ′ , Λ(l) ∧ Λ(l ′ ) ⇔ false.

• Determinism on events: ∀E,

E ′ ⊆ E with E ∩ E ′ ̸ = ∅, ∀l 0 , l, l ′ ∈ L, if l 0 E,γ,U ----→ l and l 0 E ′ ,γ ′ ,U ′ -----→ l ′ with (l, E, γ, U ) ̸ = (l ′ , E ′ , γ ′ , U ′ ), then either Λ(l) ∧ Λ(l ′ ) ⇔ false or γ ∧ γ ′ ⇔ false.
In other terms, if we are in the location l 0 and an event E 0 ∈ E ∩ E ′ occurs, then at most one transition can be fired.

• Determinism on asynchronous transitions ♯:

∀l 0 , l, l ′ ∈ L, if l 0 ♯,γ,U ---→ l and l 0 ♯,γ ′ ,U ′ ----→ l ′ with (l, γ, U ) ̸ = (l ′ , γ ′ , U ′ ), then either Λ(l) ∧ Λ(l ′ ) ⇔ false or γ ∧ γ ′ ⇔ false.
In other terms, if we are in the location l 0 , at most one asynchronous transition can be fired.

• No ♯-labelled loops: For all l 0 ---------→, there exists i ≤ p so that E i ̸ = ♯. It means the automaton cannot loop infinitely with asynchronous transitions.

E 0 ,γ 0 ,U 0 -----→ l 1 E 1 ,γ 1 ,U 1 -----→ . . . E p-1 ,γ p-1 ,U p-1 -

Hybrid Automata Stochastic Logic

Remark 4.4.2 (Events for Chemical Reaction Network)

In the case of Chemical Reaction Network, the set of events E is the set of reactions {R 1 , . . . , R M }.

Synchronised simulation

An LHA A is synchronised with the simulations of a GSPN model M. Each path is infinitely simulated when either an accepting state of A is reached, or the synchronisation can no longer proceeds (hence the path is rejected). The transitions of A synchronises with the transitions of the trajectory σ being sampled. An LHA admits two kinds of transitions.

• Synchronising transitions (associated with a subset of E, with ALL = E), which may be traversed when an event (in E) is observed.

• Asynchronous transitions (denoted by ♯), which are traversed autonomously.

Asynchronous transitions have priority over synchronised ones. Chapter 4. Verification of Continuous-Time Markov Chains

l 0 ṫ:1 l 1 ṫ:1 l 2 ṫ:1 ♯,true, {n=x I } { R 1 , R 2 } , t r u e , { n = x I } ♯ , n ≥ 8 , ∅

77

• First, the synchronised state is initialised: (m 0 , l 0 , (0.0, 0.0), 0.0).

• From l 0 , there exists only one edge l 0 → l 1 . This edge is asynchronous, so it does not need a reaction occurrence to be fired. The transition predicate is true, so l 0 → l 1 is fired, and the variable n is updated. The state becomes (m 0 , l 1 , (5.0, 0.0), 0.0).

• One asynchronous edge arises from l 1 , but the transition predicate γ = n ≥ 8 is not satisfied, the edge is not fired. The other edge is synchronous, so we have to wait that a reaction occurs.

• Time goes on until a duration of 0.2. In the location l 1 , ṫ = 1 gives that t has a flow of 1 in l 1 . More precisely, F low(l 1 ) = (0.0, 1.0) which means the derivative of t with respect to time is 1, whereas the derivative of n is zero. Thus t = 0.2.

• Then, R 1 occurs. The marking evolves to m 1 = (94, 6, 0). One edge l 1

{R 1 ,R 2 },true,{n=x I } -------------→ l 1 is synchronous. As R 1 ∈ {R 1 , R 2 } and γ = true is verified, this edge is fired.
The synchronised state becomes ((94, 6, 0), l 1 , (6.0, 0.2), 0.2).

• The transition predicate γ = n ≥ 8 is still not verified, so l 1 → l 2 is still not fired.

• The last two steps are repeated until we arrive at the time t = 0.2 + 0.1 + 0.08 + 0.02 + 0.12 = 0.52. Then R 1 occurs, and the synchronised state is updated:

((91, 8, 1), l 1 , (8.0, 0.52), 0.52). This time, the transition predicate γ = n ≥ 8 is fulfiled, so l 1 → l 2 is fired. As l 2 is a final state, the simulation ends with state ((91, 8, 1), l 2 , (8.0, 0.52), 0.52). The path is accepted by the LHA.

Note that with the path σ = (95, 5, 0) 0.12 --→ (95, 4, 1)

0.13 --→ (95, 3, 2) 0.11 --→ (95, 2, 3) 0.021 ---→ (95, 1, 4) 0.38
--→ (95, 0, 5), the LHA would have rejected σ.

HASL Expressions

An HASL formula is described by a couple (Z, A), where A is an LHA and Z is an expression built on top of the grammar in Definition 4.4.3. It leads to the evaluation of an arithmetic expression based on synchronised simulations. Such an arithmetic expression is formally described by the HASL expression Z.

Hybrid Automata Stochastic Logic Definition 4.4.3 (HASL Expression)

An HASL expression Z is defined by the following grammar

Z ::= c | P ROB | AV G(Y ) | Z + Z |Z • Z Y ::= c | Y + Y | Y • Y | Y /Y | Last(y) | M in(

Cosmos Statistical Model Checker

Cosmos [START_REF] Ballarini | Applications of an expressive statistical model checking approach to the analysis of genetic circuits[END_REF] is a statistical model checker for HASL specifications developed in C++. It takes as inputs a Generalised Stochastic Petri Net (GSPN) and an HASL formula (i.e. an LHA and an HASL expression).

This tool implements the simulation of any Discrete-Event Stochastic Processes described by a GSPN. This formalism allows efficient path simulation of the model.

Each simulation is synchronised with the LHA, and evaluation of the HASL expression is done during the simulation. For more details, we refer the reader to (Djafri, 2012, Chapter 6) and (Barbot, 2014, Chapter 6).

Script 3 (An example of a Cosmos run)

An example of a Cosmos script is available in code/chap3/example_cosmos_sir.sh with the following HASL formula that computes the average time the SIR model reaches eight infected people.

Z : AV G(Last(t))

l 0 ṫ:1 l 1 ṫ:1 l 2 ṫ:1 ♯,true, {n=x I } { R 1 , R 2 } , t r u e , { n = x I } ♯ , n ≥ 8 , ∅

Summary

In this chapter, we described:

• Temporal logics to formally describe specifications over a CTMC.

• Model checking of CTMCs with a particular focus on two problems (estimation and threshold problem) and the related statistical methods.

• Model checking on collections of CTMCs called model checking of parametric

Continuous-Time Markov Chain and the related statistical methods.

• The HASL formalism designed for discrete event stochastic processes such as CTMCs.

Summary

At this point of the thesis, we have defined the main concepts to present our contributions in statistical inference and verification of Continuous-Time Markov Chains, based on ABC methods with Linear Hybrid Automata.

Chapter 5

Automaton-ABC for the statistical inference of CTMCs

The difficulty of statistical inference for complex models relies on the lack of a closed form of the likelihood. Continuous-Time Markov Chains form a class of stochastic models that suffers from this difficulty: in most cases, the likelihood of observations is intractable.

Thus, ABC likelihood-free methods (Section 3.3) seem adequate for the Bayesian estimation of such models. In this chapter, our goal is to define the statistical framework of CTMCs and apply ABC methods using the synchronised simulation included in the HASL formalism (Section 4.4). We address several tasks:

• In Section 5.1, we describe the observation models of CTMCs. We consider two types of observations: when trajectories are observed continuously, and when trajectories are observed discretely in time.

• In Section 5.2, we detail a new algorithm called automaton-ABC, which combines classical ABC methods with synchronised simulations (defined in Section 4.4.2).

• Section 5.3 and Section 5.4 are applications of the automaton-ABC algorithm for two different tasks: take into account oscillatory trends of models and Bayesian inference.

Observation model and likelihood

Let us consider observations y exp . A statistical model induced by a parametric CTMC (Definition 4.3.1) is a dominated model [START_REF] Billingsley | Statistical Inference for Markov Processes[END_REF]) (Proposition 5.1.1). Let S a CTMC. P S admits a density p S :

Observation model and likelihood

p S :P ath(M) → R ≥0 σ = s 0 t 0 -→ . . . t k-1 ---→ s k → P (s 0 ) k i=1 Q(s i-1 , s i )e -E(s i-1 )t i-1
Such properties are also helpful in the Bayesian framework (Proposition 5.1.2).

Proposition 5.1.2 (Conditional density of a parametric CTMC)

Let (S θ ) θ∈Θ be a parametric CTMC, P ϑ a prior distribution.

• (θ, A) → P S θ (A) is a transition kernel.

• P S |ϑ=θ admits a continuous density σ → p(σ|θ) = p S θ (σ).

Before any mention of learning algorithm from data, we need to formulate the nature of observations y exp mathematically in the case of Continuous-Time Markov

Chains. We consider two types of observations.

• Event-discrete observations. We observe a path σ ∈ P ath(M) by observing the transitions of the CTMC.

• Time-discrete observations. At fixed times t 1 , . . . , t N , we observe the state of the CTMC. It corresponds to the classical case of observation of a state-space dynamical model.

Event-discrete observations

Let S θ = (S t θ ) t∈R ≥0 be a CTMC defined by a transition rate matrix Q θ , transition probability matrix P θ , and exit rate E θ . We consider the paths σ (1) , . . . , σ (N ) as observations.

Definition 5.1.1 (Likelihood function of a parametric CTMC)

Let σ (1) , . . . , σ (N ) observation paths. The likelihood function is defined as:

l(θ; σ (1) , . . . , σ (N ) ) = N k=1 p(σ (i) |θ) = N i=1 P (σ (i) [0]) k (i) j=1 Q θ (σ (i) [j-1], σ (i) [j])e -E θ (σ (i) [j-1]δ(σ (i) ,j-1))
where k (i) is the number of states of σ (i) .

These observations correspond to the case where we fully observe the transitions of the CTMC model. This case seems unrealistic for most biological phenomena, but it is more reasonable for queueing network models or fault tolerance models because computers can log any occurred event. An analytical form of the Maximum Likelihood

Estimator (Definition 3.1.3) can be derived (Bladt and Nielsen, 2017, Chapter 12).

Sometimes, we can only observe a subset of the state space. Thus, we introduce an observation function γ defined over the set of trajectories by:

γ : P ath(M) → E σ → γ(σ) = g(s 0 ) t 0 -→ . . . t k -→ g(s k )
where g is the projection of an element s ∈ S over the set of observed states. Thus, γ(σ) is a trajectory with a reduced state space.

A frequent case is when only a few variables of the system are observed. In the Chemical Reaction Network framework, it means only a few species are observed. In this case, we define I g = {i 1 , . . . , i d obs } ⊆ {1, . . . , d} as the set of observed species and the observation function g is:

g : S → N d obs s → (s[i 1 ], . . . , s[i d obs ])
For example, with the SIR model 2.15, we usually only observe the number of infected people, so the observation function over the state space g is defined as

g : N 3 → N (X S , X I , X R ) → X I .

Time-discrete observations: state-space model

In most applications of continuous phenomena, observations are considered discrete in time. In this case, the standard framework used for the observation model is the state-space model [START_REF] Cappé | Inference in Hidden Markov Models[END_REF].

Let S = (S t ) t∈R ≥0 a CTMC. Let t 1 , . . . , t K a collection of observation times.

We suppose Z = (S k ) 1≤k≤K is unobservable, and we consider Y = (Y k ) 1≤k≤K the corresponding observable states. The state-space model is the Hidden Markov Model (Y, Z) defined by the scheme, for k ∈ {1, . . . , K}:

Y k ∼ g(.|S t k , θ y ) Z k = S t k ∼ p(.|S v , θ) with t 0 ≤ v ≤ t k .
(Y k ) k∈ N are the observations of some hidden/latent continuous trajectory. g is the observation function, it models the link between the observed state and the latent state. In this case, the i.i.d random observations are named (Y (1) , . . . , Y (N ) ), which a realisation is (y

(1) 1:K , . . . , y (N ) 1:K ).
Due to the Markov property of CTMCs, the discrete scheme is rewritten as:

Y k ∼ g(.|S t k , θ y ) S t k ∼ p(.|S t k-1 , θ) .
Let us compute the likelihood for one observation y 1:K ∈ S K :

l(θ; y 1:K ) = p Y (y 1:K ) = z 0:K ∈S K+1 p Y,Z (y 1:K )dz 0:K = z 0:K ∈S K+1 p Y |Z (y 1:K | z 0:K )p Z (z 0:K )dz 0:K = z 0:K ∈S K+1 p(z 0 ) N k=1 g(y k |z k )p(z k |z k-1 )dz 0:K For i.i.

d observations y

(1) 1:K , . . . , y (N ) 1:K , the likelihood is:

l(θ; y (1) 1:K , . . . , y (N ) 1:K ) = N i=1 z 0:K ∈S K+1 p(z 0 ) N k=1 g(y (i) k |z k )p(z k |z k-1 )
1:K ) = N i=1 p(y 0 ) N k=1 p(y (i) k |y (i) k-1 ).
Even in this simple case, the likelihood is difficult to compute. Indeed, the transition function of a CTMC p(y k |y k-1 ), i.e. the probability of being in a state at time t k given the state at time t k-1 , is challenging to compute. The observations times are not coincident with jumps/transitions of the system, so we only know

p(y k |y k-1 ) = P(S t k = y k |S t k-1 = y k-1 ) = δ y k-1 (y k )+Q(y k-1 , y k )(t k -t k-1 )+o(t k -t k-1 ).
Note that given observations y k , t k and a CTMC, methods that reconstruct the most probable trajectory σ exist, such as Viterbi algorithm or more evolved methods [START_REF] Perkins | Maximum likelihood trajectories for continuous-time Markov chains[END_REF].

Statistical inference of Continuous-Time Markov Chains based on time-discrete

observations has been much studied in the literature [START_REF] Craciun | Statistical model for biochemical network inference[END_REF][START_REF] Schnoerr | Approximation and inference methods for stochastic biochemical kinetics -A tutorial review[END_REF][START_REF] Loskot | Comprehensive review of models and methods for inferences in bio-chemical reaction networks[END_REF]. Consideration of time-discrete observations for CTMCs has been considered in [START_REF] Bladt | Statistical inference for discretely observed Markov jump processes[END_REF]. They proved the existence and uniqueness of the maximum likelihood estimator under certain conditions, and the possibility of applying Expectation-Maximisation (EM) and MCMC procedures. More evolved EM algorithms suited to CTMCs exist. For example, stochastic approximation EM algorithm that involves an ABC procedure is detailed in [START_REF] Picchini | Coupling stochastic EM and approximate Bayesian computation for parameter inference in state-space models[END_REF]. Some of these statistical methods have been implemented in R [START_REF] Pfeuffer | Ctmcd: An R package for estimating the parameters of a continuous-time Markov chain from discrete-time data[END_REF].

Particularly, Bayesian inference for CTMCs has been considered [START_REF] Gómez-Corral | Bayesian Inference of Markov Processes[END_REF]. [START_REF] Warne | Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art[END_REF] reviews the main Bayesian techniques for CTMCs with a particular focus on ABC related methods. They are illustrated by numerical examples.

In the thesis [START_REF] Alharbi | Bayesian Inference for Continuous Time Markov Chains[END_REF], Bayesian inference is considered on CTMCs, focusing on the two main likelihood-free methods of the Bayesian statistical literature: ABC and Particle Marginal Metropolis-Hasting algorithm [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF]. Note that to compare the likelihood-free posteriors, in (Alharbi, 2018, Chapter 6) More recently, ABC methods have been coupled with Multi-Level Monte Carlo for CTMC inference. The main idea of Multi-Level Monte Carlo is to estimate (the expectation of) a probability distribution with approximate samples. These samples have different levels of accuracy l ∈ {1, . . . , L}. The Multi-Level Monte Carlo estimator reduces the variance compared to the classical Monte Carlo estimator. In [START_REF] Jasra | Multilevel Monte Carlo in approximate Bayesian computation[END_REF], the approximation is related to an ABC tolerance ϵ l , l ∈ {1, . . . , L}, i.e. the approximate samples are drawn from π ϵ l ABC , {1, . . . , L}. In [START_REF] Lester | Multi-level Approximate Bayesian Computation[END_REF], the approximation is related to the accuracy of the CTMC simulation. Instead of exact CTMC simulation (Stochastic Simulation Algorithm, see Section 2.3.4), simulation is performed by the Tauleap Algorithm 3 with a time step τ l , l ∈ {1, . . . , L}.

ABC methods have proven efficiency in the context of CTMCs. In this chapter, we will consider ABC procedures with a verification point of view for Bayesian inference.

Approximate Bayesian Computation for event-discrete observations

As we saw in Section 3.3 and above in this chapter, Approximate Bayesian Computation methods are suitable for Bayesian inference when the likelihood is complex.

Hence, these methods are helpful for Continuous-Time Markov Chains. In this section, we detail several configurations of ABC methods in order to show how they are usable, and apply them to several models.

Distance over paths of CTMC

Here, we consider event-discrete observations (Section 5.1.1) y exp = (υ (1) , . . . , υ (n) )

where

υ (i) = g(s 0,(i) ) t 0,(i) ---→ . . . t k (i) -1,(i) ------→ g(s k (i) ,(i) ), I g ⊆ {1, .
. . , d} is the set of observed species, and g ∈ {0, 1} d the observation function. By denoting g(s j,(i) ) = y (i) j , one can represent the i-th observation with (y

(i) j , t (i) j ) j∈0:k (i) .
The observation paths are defined over [0, T end ], T end ∈ R ≥0 , which means j t j,(i) = T end .

To execute the ABC algorithm, we have to define a distance function between observations and simulations. A classical choice is the Euclidean distance:

d ||.|| 2 : R k × R k → R ≥0 (y (1) , y (2) ) → j (y (1) j -y (2) j ) 2
(5.1) However, the dimension of our observations is not fixed. To use the Euclidean distance, we have to define a discrete timeline t 1 , . . . , t N over [0, T end ] and collect the values of the paths γ(σ (i) ) over this timeline. This operation is equivalent to consider the scheme of time-discrete observations (Section 5.1.2), which lose information.

As the observations are step functions on a bounded interval [0, T end ], one can use L 1 and L 2 distances for 1-dimensional trajectories:

d L p (σ 1 , σ 2 ) = ( [0,T end ] |σ 1 @t -σ 2 @t| p dt) 1 p , p ∈ {1, 2} (5.2)
which guarantees the full use of the information contained in the observations.

As our observations are step functions, the integral is easy to compute. Fig- ure 5.1 illustrates the computation for two one-dimensional paths σ 1 and σ 2 . The distance between the paths is the sum of the green areas. Hence, the computational cost of the distance is proportional to the number of events for both trajectories.

Observation model and likelihood

Examples of ABC-SMC inference on parametric CTMC with different observation schemes and distances

In this section, we present a few examples of Bayesian inference with ABC-SMC (Algorithm 14). We compare the effect of three observation configurations based on the same simulated dataset. We see the evolution of the resulting ABC posteriors around the true parameter estimated with a fixed number of simulations.

We consider a simulated dataset of five observations y exp = (υ (1) , . . . , υ (5) ) and the distance:

ρ : E 5 × E 5 → R ≥0 (y, y exp ) → min 1≤i,j≤5 d(y[i], y exp [j] υ (j) ) (5.3)
where d can either be the integral L 2 distance (5.2) or the Euclidean distance (5.1).

For each model, two sets of observation times are used in the case of Euclidean distance. Observations are collected according to these observation times.

Resulting posterior distributions should be comparable, so we fix a maximal number of simulations N end ∈ {10 6 , 10 7 }. N end is the simulation budget. When the number of simulations reaches N end , we let the current step of the ABC-SMC algorithm finish and end the algorithm's run. In practice, it results in a similar number of simulations. In this application, the adaptive tolerance schedule is set with α = 0.75 (Section 3.3.2).

Enzymatic reaction model

We consider the enzymatic reaction (ER) system. A substrate species S is converted into a product P through the mediation of an enzyme E. The processes are modelled by the CRN (5.4), parametrised by the parameter vec-

tor θ = (k 1 , k 2 , k 3 ). k 1 , k 2 , k 3 are the kinetic rate constants of the reactions R 1 , R 2 , R 3 .
The initial state is (E 0 , S 0 , ES 0 , P 0 ) = (100, 100, 0, 0).

R 1 : E + S k 1 -→ ES R 2 : ES k 2 -→ E + S R 3 : ES k 3 -→ E + P (5.4)
We simulate the dataset (five observations) with θ true = (1, 1, 5) and T end = 0.075. Parameter inference focuses on the parameter k 3 (the two others are supposed known) with the prior distribution k 3 ∼ U (0, 100).

Three observation settings are considered with an end time of T end = 0.075:

Chapter 5. Automaton-ABC for the statistical inference of CTMCs

89

• 5 event-discrete observations y exp = (υ (1) , . . . , υ (5) ), with L 2 distance.

• 5 time-discrete observations with 10 points y (1) , . . . , y (5) with y (i) = (y

(i) 1 , . . . , y (i)
10 ) collected from (υ (1) , . . . , υ (5) ), with Euclidean distance.

• 5 time-discrete observations with 75 points y (1) , . . . , y (5) with y (i) = (y

(i) 1 , . . . , y (i)

75

) collected from (υ (1) , . . . , υ (5) ), with Euclidean distance.

One of the trajectories from the simulated dataset is shown in Figure 5.2 with the two time-discrete observations settings.

SIR model

We consider the SIR model already described in 2.15 with initial state s 0 = (95, 5, 0). We simulate the dataset with θ true = (1.2E-3, 5.0E-2) and T end = 150.0. Both parameters are infered with prior distributions k i ∼ U (5E-5, 3E-3) and

k r ∼ U (5E-3, 0.2).
The three observations settings are:

• 5 event-discrete observations y exp = (υ (1) , . . . , υ (5) ), with L 2 distance.

• 5 time-discrete observations with 10 points y (1) , . . . , y (5) with y (i) = (y

(i) 1 , . . . , y (i)
10 ) collected from (υ (1) , . . . , υ (5) ), with Euclidean distance.

• 5 time-discrete observations with 150 points y (1) , . . . , y (5) with y (i) = (y

(i) 1 , . . . , y (i) 150 )
collected from (υ (1) , . . . , υ (5) ), with Euclidean distance.

One of the trajectories from the simulated dataset is shown in Figure 5.4 with the two time-discrete observations settings.

d ||.|| 2 , 10 points d ||.|| 2 , 150 points d L 2 Mean k i 1.27E-3 1.32E-3 1.31E-3 Std k i 3.31E-4 3.14E-4 3.14E-4 Mean k r 4.70E-2 4.75E-2 4.72E-2
Std k r 7.54E-3 6.35E-3 6.35E-3

(k i , k r ) = (1.2E-3, 5E-2).

Discussion

The most time-consuming experiment (SIR model with N end = 10 7 , 150-points observations) has not exceeded 35 minutes with eight jobs. These two experiments showed that adding information in the observations does not constantly improve the inference in likelihood-free ABC methods. Indeed, the little information contained in the addition of more observations is counterbalanced by the curse of dimensionality. Adding more points leads to the increase of the observation space 5.2. Automaton-ABC: ABC procedures with synchronised simulation dimension. Thus, an important task is to compute informative statistics (Sisson, Fan, and Beaumont, 2019, Chapter 5) over the observations. This task can be considered more important than collecting more observations in some cases. An original approach with parametric CTMC oscillators is presented in Section 5.3.

Script 4 (Examples of ABC inference on parametric CTMC)

The scripts of these ABC inference examples can be found in the git repository of the thesis at code/chap4/inference.

Automaton-ABC: ABC procedures with synchronised simulation

In Section 4.4.2 of Chapter 4, we have defined synchronised simulations of CTMCs with a Linear Hybrid Automaton.

Evaluating an HASL trajectory expression Y with the synchronised simulation of a CTMC can be seen as an efficient method to compute complex summary statistics based on properties. Indeed, the automaton variables are updated within the simulation, which can be seen as statistics of the simulated trajectory. In this context, HASL formalism offers two main advantages:

• It offers a well-suited and visual tool for the definition of properties over CTMC trajectories.

• A synchronised simulation does not need the storage of trajectory values which saves memory allocation time.

In the following, we formulate an adaptation of ABC procedures (Section 3.

(i) ) 1≤i≤N drawn from π ϵ ABC for i = 1 : N do repeat θ ′ ∼ π d ′ ∼ (Y, A) × M θ ′ until d ′ ≤ ϵ θ (i) ← θ ′ end for
Algorithm 16 is the Sequential Monte Carlo version of Algorithm 15, with the adaptive tolerance schedule of Section 3.3.2.

Algorithm 16 General automaton-ABC Sequential Monte Carlo Algorithm

Require: (M θ ) θ∈Θ parametric CTMC, π prior, N : number of particles, ϵ: tolerance level, α: quantile level, Y HASL trajectory expression, A LHA, K: kernel proposal density Ensure:

(ω (i) , θ (i) ) 1≤i≤N weighted samples drawn from π ϵ ABC θ (i) 0 ∼ π, i ∈ 1, . . . , N d i ∼ (Y, A) × M θ (i) 0 , 0 ∈ 1, . . . , N ϵ current ← quantile(α, (d i ) 1≤i≤N ) (ω (i) 0 ) 1≤i≤N ← 1 N m ← 1 while ϵ current > ϵ do for i = 1 : N do repeat Take θ ′ from (θ (j)
m-1 ) 1≤j≤N with probabilities (ω

(j) m-1 ) 1≤j≤N θ (i) m ∼ K(.|θ ′ ) d i ∼ (Y, A) × M θ (i) m until d i ≤ ϵ current ω (i) m ← π ( θ (i) m ) N Σ i ′ =1 ω (i ′ ) m-1 K(θ (i) m |θ (i ′ ) m-1 )
end for Normalise (ω

(i) m ) 1≤i≤N ϵ current ← quantile(α, (d i ) 1≤i≤N ) m ← m + 1 end while return (ω (i) m , θ (i) m ) 1≤i≤N
The fact that y exp is not mentioned does not mean that these procedures cannot be related to statistical inference. Indeed, the constants of the automaton A may 5.3. Oscillatory trends of genetic networks contain the observations, which will be the case in Section 5.4. In the rest of the chapter, as well as the next chapter, automaton-ABC procedures are used to tackle several problems of statistical inference and statistical verification.

Oscillatory trends of genetic networks

In computational biology, quantitative analysis of the biological system's dynamics is challenging even for simple models [START_REF] Baldan | Petri nets for modelling metabolic pathways[END_REF].

This section focuses on a specific behaviour of dynamical systems: oscillatory trends of gene regulatory networks. Genetics oscillators have grown much interest in the domain of synthetic biology over the past two decades [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF][START_REF] O'brien | Modeling synthetic gene oscillators[END_REF][START_REF] Li | Systems and synthetic biology approaches in understanding biological oscillators[END_REF]. They allow modelling many rhythmed physiological phenomena such as circadian clock [START_REF] Leloup | Toward a detailed computational model for the mammalian circadian clock[END_REF], heartbeat or cell cycle.

Characterisation of oscillations in model checking CTMC has been considered in the literature [START_REF] Ballarini | Analysing Biochemical Oscillation through Probabilistic Model Checking[END_REF] (Spieler, 2014, Chapter 7). In [START_REF] Ballarini | Applications of an expressive statistical model checking approach to the analysis of genetic circuits[END_REF], they use the expressiveness of HASL for the characterisation of noisy oscillations. In this work, we define a Linear Hybrid Automaton A per based on [START_REF] Ballarini | Applications of an expressive statistical model checking approach to the analysis of genetic circuits[END_REF] that measures a CTMC trajectory time period and computes a distance to oscillations. We explore the parameter space of parametric CTMCs that show oscillatory behaviours with an indicated time period, based on the LHA A per . The main idea is to use this automaton that measures the mean time period of a trajectory within automaton-ABC (Section 5.2) to find the subset of the parameter space that produces trajectories with an indicated time period.

Period automaton

A per l 0 begin l ′ 0 ṫ = 1 low ṫ = 1 ṫp = 1 mid ṫ = 1 ṫp = 1 high ṫ = 1 ṫp = 1 final x O > L x O ≤ L L < x O < H x O ≥ H ALL, true, ∅ ALL, true, ∅ ALL, n<N, ∅ ALL, n<N, ∅ ALL, n<N, ∅ ♯, t ≥ initT, ∅ A L L , t r u e , { t p : = 0 , n = - 1 } ALL, true, {tp:=0, n= -1} ALL, n<N, ∅ ALL, n < N, {top := 1} ALL, n<N, ∅ ALL, (n<N ∧ top = 0), ∅ ALL, (n= -1 ∧ top = 1), {n++, tp:=0, top:=0} ALL, (0 ≤ n ≤ 1 ∧ top=1), {n++, top:=0, m( tp, tp, n), tp:=0} ALL, (2 ≤ n < N ∧ top=1), {n++, top:=0, tp=m( tp, tp, n), ft p =v(ft p , tp, tp, n), tp:=0} ♯, n=N, d:=dist( tp, ft p , tp (obs)
)

♯, n=N, d:=dist( tp, ft p , tp (obs) ) ♯ , n = N , d : = d i s t ( t p , f t p , t p ( o b s ) ) Figure 5
.6: Period automaton Aper.

We consider the LHA A per in Figure 5.6 to measure oscillatory trends of trajectories.

This automaton is designed to compute two measures for the trajectory of a species O: the period duration mean tp and the period variance f tp .

It is mainly formed of three locations: low, mid and high, plus other locations related to the initialisation and the end of the simulation. These three locations outline three levels of the species population O. Being in the location low means that the species population is below a user-defined level L, i.e. x O ≤ L, whereas entering mid (resp. high) means L < x O < H (resp. x O ≥ H). L and H are hyperparameters of the automaton and have to be defined during the creation of the automaton. Moreover:

• n is the number of detected periods.

• top is a boolean that indicates if the automaton has entered high location.

• tp is the period mean, and f tp is the period variance. 96

Oscillatory trends of genetic networks

• t p is the current computed time period. It equals zero when tp and f tp are updated after the detection of a period, i.e. when mid → low occurs with top ⇔ true.

• L, H and N per are hyperparameters of the automaton. L, H defines the partition of the population species. N per is the number of periods to be detected. After the detection of N per periods, the automaton ends.

• d is the distance from a period mean reference called tp (obs) .

A period duration is computed each time a trajectory contains a sub-trajectory that corresponds to a Aper trace of the form low

• (low|mid) * • high • (high|mid) * • mid • low:
• We start to count a period when the automaton enters low.

• The trajectory oscillates between low and mid.

• The automaton enters high at least once.

• The period duration finishes when mid → low is observed with top ⇔ true. • The grey part corresponds to the initial locations l 0 and l ′ 0 . The automaton is in one of these locations until it reaches low, i.e. x O ≤ L.

• When the automaton reaches low, n = -1. The automaton waits to detect a period before computing the measures.

• The trajectory continues until a period is detected, n is updated to 0.

• The trajectory continues until another period is detected. n is updated to one, and the period mean is updated to t p .

• For the next period, the period mean and variance period are updated according to Equations 5.5.

• The trajectory continues until n reaches N per = 5.

m( tp , t p , n) = (n -1) tp + t p n v(f tp , tp , t p , n) = n -2 n -1 • f tp + ( tp -t p ) 2 n (5.5)
At the end of the simulation (i.e N per periods are detected), the variable d is updated. The automaton computes a distance from a duration period mean.

dist( tp , f tp , tp (obs) ) = min( tp -tp (obs) tp (obs)
, f tp tp (obs) )

The variation coefficient helps discard trajectories with, for example, periods of duration 0.5 tp (obs) and 2 tp (obs) but a resulting period mean of tp (obs) . Thus, if d ≤ 0.1, the relative error from the period mean reference is less than 10%, as well as the variation coefficient.

Applications of the automaton-ABC algorithm with Aper

In the following, we use the automaton-ABC algorithm with the period automaton Aper.

Doping 3-way oscillator

Model. First, we consider a toy model for oscillatons called doping 3-way oscillator, from which a probabilistic model checking description has been detailed in [START_REF] Ballarini | Analysing Biochemical Oscillation through Probabilistic Model Checking[END_REF]. It is described by the following six reactions with massaction law:

Oscillatory trends of genetic networks

A + B r A -→ 2B B + C r B -→ 2C C + A r C -→ 2A (5.6) D A + C r C -→ A D B + A r A -→ D B + B D C + B r B -→ D C + C (5.7)
The number of species is constant in time. Equations 5.6 creates oscillations of species by creating positive loop feedback. If the first reaction occurs a few times, the number of species B increases, which makes the second reaction (that induces consumption of B) more probable, and so on. Equations 5.7, called doping reactions, permit sustained oscillations by avoiding that one species vanishes.

Results. We apply the automaton-ABC algorithm with the synchronised simulation (Last( tp ), Aper). We want to find the subset of parameters for which oscillations have a period duration mean of 0.02.

We made two experiments. The first is one-dimensional, the setting is:

• N per = 4 periods, L = 300 and H = 360

• tp (obs) = 0.01

• r A ∼ U (0, 10) and other parameters are fixed (r B = r C = 1.0). 333,333,333,10,10,10) • N = 1000 particles

• s 0 = (A 0 , B 0 , C 0 , (D A ) 0 , (D B ) 0 , (D C ) 0 ) = (
• Tolerance of 20% (ϵ = 0.2) we have reduced the parameter space to a subset where it is probable to obtain trajectories with a period mean of 0.01 (relatively to a tolerance of 20%). Second, the posterior has only one mode, which is quite intuitive. As we have fixed r B and r C , one could expect that the period duration mean is directly linked to the kinetics of reaction 1, which is only parametrised by r A . The second experiment considers variations over the whole parameter space. The setting is:

• N per = 4 periods, L = 300 and H = 360

• tp (obs) = 0.01

• r A , r B , r C ∼ U (0, 10) 333,333,333,10,10,10) • N = 1000 particles

• s 0 = (A 0 , B 0 , C 0 , (D A ) 0 , (D B ) 0 , (D C ) 0 ) = (
• Tolerance of 20% (ϵ = 0.2) Figure 5.9 shows the correlation plot of the resulting automaton-ABC posterior.

The parallel run (120 jobs) performed 533003 simulations and lasted 268 seconds.

We can see that the different parameters have structured correlations: the three histograms have parabolic shapes. Also, one can notice that for each histogram, the area near the point (1, 1) is a high probability area. This is consistent with the above one-dimensional experiment.

Repressilator model

Model. Repressilator models are synthetic genetic networks that describe a cycle of three genes with a negative feedback loop [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF]. In this part, we consider one of the first repressilator model developed by [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. This model includes three genes lacI, tetR and cI (called G i , i ∈ {1, 2, 3} in the equations) and their corresponding proteins LacI, TetR and CI (called Equations 5.8, 5.9 and 5.10 below constitute a possible reduced description of the processes involved in the repressilator dynamics.

P i , i ∈ {1, 2, 3}). mRN A i is
tr1 : G 1 r 1 -→ G 1 + mRN A 1 tr2 : G 2 r 2 -→ G 2 + mRN A 2 tr3 : G 3 r 3 -→ G 3 + mRN A 3 (5.8) trl1 : mRN A 1 β -→ mRN A 1 + P 1 trl2 : mRN A 2 β -→ mRN A 2 + P 2 trl3 : mRN A 3 β -→ mRN A 3 + P 3 (5.9) degr1 : mRN A 1 1 -→ ∅ degr2 : mRN A 2 1 -→ ∅ degr3 : mRN A 3 1 -→ ∅ degr4 : P 1 1 -→ ∅ degr5 : P 2 1 -→ ∅ degr6 : P 3 1 -→ ∅ (5.10)
The kinetic rates are:

• r 1 = α

1+[P 3 ] n + α 0 , r 2 = α 1+[P 1 ] n + α 0 and r 3 = α 1+[P 2 ] n + α 0 for the transcription reactions,
• mass-action law for the translation reactions (i.e. β • [mRN A i ]),

• 1 for the degradation reactions (independent of the concentrations of P i ).

These reaction channels do not follow mass-action law. Indeed they use a reparametrisation of the model detailed in (Elowitz, Leibler, and Leibier, 2000, Box 1), which is convenient to study the dynamics of the related Reaction Rate Equations of the CRN (continuous ODE system):

• β is the ratio of the protein decay rate over the mRNA decay rate. Equations 5.9 are thus parametrised by β, whereas degradation Equations 5.10 are not parametrised.

• n is the Hill coefficient: it measures the interaction affinity between species, i.e. how much the ligands will bind to create macromolecules.

• α is the parameter related to transcription growth.

• α 0 is the parameter related to a minimum level of transcription growth. Indeed

r i → [P i+1 ]→+∞
α 0 , which means α 0 models a saturation of the repressors' effect.

With Equations 5.8 5.9 5.10, we neglect several aspects of molecular interactions:

• Binding and unbinding of RNA polymerase with promoters, which is a necessary phenomenon for any transcription, is implicit in Equations 5.8.

Oscillatory trends of genetic networks

• The observation of one protein species is initially measured in [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] by a green fluorescent protein (GFP), which has the role of a reporter for the TetR protein. Our model supposes that we directly observe the protein (or the link between GFP and TetR is known and invertible).

Simulations. To better understand the system's dynamics, we present a few simulations. Default parameter values are (α, β, n, α 0 ) = (200.0, 2.0, 2.0, 0.0).

Oscillatory trends of genetic networks

Results. We apply the automaton-ABC algorithm with the synchronised simulation (Last(d), Aper). We want to find the subset of parameters for which oscillations have a period duration mean of 20.0.

We did two experiments. The first is three-dimensional, the setting is:

• N per = 4 periods, L = 50 and H = 200

• tp (obs) = 20.0

• α ∼ U (50, 5000), β ∼ U (0.1, 5.0), n ∼ U (0.5, 5.0) (α 0 = 0.0).

• s 0 = ((mRN A 1 ) 0 , (mRN A 2 ) 0 , (mRN A 3 ) 0 , (P 1 ) 0 , (P 2 ) 0 , (P 3 ) 0 ) = (0, 0, 0, 5, 0, 15)

• N = 1000 particles

• Tolerance of 10% (ϵ = 0.1) The second experiment is four-dimensional. The setting is:

• N per = 4 periods, L = 50 and H = 200

Chapter 5. Automaton-ABC for the statistical inference of CTMCs

105

• tp (obs) = 20.0

• α ∼ U (50, 5000), β ∼ U (0.1, 5.0), n ∼ U (0.5, 5.0), α 0 ∼ U (0.0, 5.0).

• s 0 = ((mRN A 1 ) 0 , (mRN A 2 ) 0 , (mRN A 3 ) 0 , (P 1 ) 0 , (P 2 ) 0 , (P 3 ) 0 ) = (0, 0, 0, 5, 0, 15)

• N = 1000 particles

• Tolerance of 10% (ϵ = 0.1)

Figure 5.16 shows the correlation plot of the resulting automaton-ABC posterior.

The parallel run (250 jobs) performed 22291 simulations and lasted 491 seconds. One can see adding a degree of freedom on α 0 has changed the correlation between α and β as well as α and n, whereas the correlation between β and n seems to have the same shape.

Script 5 (Oscillatory trends of parametric CTMC)

The scripts about oscillations of the doping 3-way oscillator and repressilator can be found in the git repository of the thesis at code/chap4/oscillators.

Accelerating the ABC procedure with HASL formalism

This section presents an automaton that aims to speed up the general ABC procedure used in Bayesian statistical inference.

In complex models, one simulation can be quite computationally expensive depending on the time window. The repressilator model is an illustrative example (cf Section 5.3.2). As ABC methods are based on a lot of model simulations, the convergence of ABC methods is directly linked to the cost of simulations.

ABC procedures simulate the model until the distance between a simulation and the observations are lower than a tolerance ϵ. However, this tolerance can be reached before the simulation ends.

Let (t i , y i ) 1≤i≤K be time-discrete observations of a one-dimensional trajectory and let σ be a path of a CTMC model. The Euclidean distance between the observations and the simulation is:

K i=1 (y i -σ@t i ) 2
If the simulation σ differs a lot from the observations, it may exist i * < K

with i * i=1 (y i -σ@t i ) 2 > ϵ.
One can deduce that, within the ABC procedure, the simulation could have been stopped at time t i * instead of simulating until t K . This section presents a Linear Hybrid Automaton called A ABC,ϵ that computes the Euclidean distance between a trajectory and observations. The synchronised simulation with A ABC,ϵ calculates the Euclidean distance to the observations on the fly of the simulation. It terminates if the distance exceeds the tolerance ϵ or time reaches t K . The main goal is to save computation time in ABC procedures by cutting off useless simulation steps when the simulation already exceeds the tolerance ϵ before t K .

Automaton

A ABC,ϵ l 0 start l 1 ṫ = 1 l 2 ♯, true, {n = x 0 , d = 0, idx = 1} ♯, t ≥ t idx , {d+ = (n -y idx ) 2 , idx+ = 1} ALL, true, {n = x O } ♯, idx ≥ Nobs + 1, {d = √ d} ♯, d > ϵ 2 , {d = +∞ } Figure 5
.17: Automaton A ABC,ϵ .

Figure 5.17 shows the A ABC,ϵ Linear Hybrid Automaton. The automaton variables are:

• n is the population of the observed species O,

• (t i , y i ) 1≤i≤K are the one-dimensional time-discrete observations,

• idx ∈ {1, . . . , K} is the index of the next observation (y idx ) to add,

• ϵ is the tolerance level of the computed distance.

For any synchronised simulation, first, the distance d is initialised to 0, and idx is initialised to 1 because the next observation to add is y 1 . It corresponds to the firing of the transition l 0 → l 1 .

While in l 1 , for any occurred reaction, n is updated (synchronous bottom self loop). If time t reaches t idx , the distance d is updated by adding (n -y idx ) 2 and idx is incremented by one (asynchronous top self edge loop). Note that while in l 1 , d is the square of the Euclidean distance.

The simulation finishes in l 2 in two ways. Either the distance has exceeded the tolerance ϵ (d > ϵ 2 ) and bottom edge l 1 → l 2 fires, or all the observations have been counted and l 1 → l 2 top edge fires.

Applications

To illustrate the potential gain of computational time for automaton-ABC with A ABC,ϵ , we consider a challenging task of inference with the repressilator model inspired from data exposed in [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF].

Accelerating the ABC procedure with HASL formalism

We considerer a simulated dataset of a 40-points observation of the population of the species P 1 from t 1 = 0.0 to t 40 = 400.0 with parameters (α, α 0 , β, n) = (200.0, 0.0, 5.0, 4.0), which is plotted in Figure 5.18. In the following, we run two parallel ABC-like procedures with 80 jobs on HPC resources from the "Mésocentre":

• Classical ABC-SMC: at each iteration, we simulate until t K and then computes the Euclidean distance from the observations.

• Automaton-ABC-SMC with (Last(d), A ABC,ϵ ).

For both ABC and automaton-ABC runs, the setting is:

• N = 500 particles,

• β ∼ U (0.5, 10.0), n ∼ U (0.5, 10.0) (α = 200.0 and α 0 = 0.0 are fixed),

• Tolerance of 50% relative error ϵ = 0.5||y exp ||. With a longer timeline of observations or a more complex model, the saved computational time could be even greater.

Script 6 (Example CTMC)

This experiment can be found in the git repository of the thesis at code/chap4/ reject_abc_automaton.

Summary

In this Chapter:

• We have detailed the statistical framework of CTMCs, the different types of observations, and distances between observations. We have discussed the available techniques for statistical inference of CTMCs.

• We have detailed automaton-ABC, a new ABC procedure based on synchronised simulations of a CTMC with a Linear Hybrid Automaton. Thus, the last statement supposes that one can measure the discrepancy between a simulation and a specification φ. We formally define a satisfiability distance of a trajectory of a CTMC from a logical property φ and use Linear Hybrid Automata to compute this distance. This distance is integrated within the ABC procedure to estimate the subset of the parameter space in which the logical property can be satisfied thanks to the obtained ABC posterior. We show that the sequential version of ABC (Algorithm 18) is well suited to using the distance from a property,

S µ × [t 1 , t 2 ],
with S µ a satisfiability region.

For example, for a bi-dimensional CTMC, the formula [[5, 8]] ×

µ 1 = x 1 ≥ 1 ∧ x 1 ≤ 2 induces the satisfiability region R 1 = [[1, 2]] × N while µ 2 = [(x 1 ≤ 3) ∨ (x 1 ≥ 5 ∧ x 1 ≤ 8)] ∧ x 2 > 4 induces the satisfiability region R 2 = ([[0, 3]] × [[5, +∞[) ∪ (
[ [5, +∞[). By a slight abuse of vocabulary, we identify µ with its satisfiability region, and we say that two formulae µ 1 , µ 2 are disjoint if their corresponding regions are.

Let µ be a propositional formula. We introduce the notion of satisfiability distance of a trajectory from a temporal formula build on top of µ over [t

t last = min{t ∈ R ≥0 , ∀t ′ ∈ [t, t 2 ]
, σ@t ′ = σ@t}. We define the distance d(σ, φ) for the following kinds of temporal formulae φ built on top of an elementary propositional formula µ:

• φ = F [t 1 ,t 2 ] µ d(σ, F [t 1 ,t 2 ] µ) =    d e ((t last , σ@t last ), S µ [t 1 ,t 2 ] ) if t last < t 1 and σ@t last ̸ ∈ S µ min t∈[t 1 ,t 2 ]
d e (σ@t, S µ ) elsewhere (6.1)

For non-elementary propositional formulae µ = µ i , we define the distance:

d(σ, F [t 1 ,t 2 ] µ i ) = min i d(σ, F [t 1 ,t 2 ] µ i )
Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs 115 where µ i are elementary formulae.

• φ = G [t 1 ,t 2 ] µ d(σ, G [t 1 ,t 2 ] µ) = t 2 t 1 d e (σ@t, S µ )dt (6.2)
Similarly, for non-elementary propositional formulae, we define the distance

d(σ, G [t 1 ,t 2 ] µ i ) = min i d(σ, G [t 1 ,t 2 ] µ i )
where µ i are elementary formulae.

• φ = µ 1 U [t 1 ,t 2 ] µ 2 d(σ, µ 1 U [t 1 ,t 2 ] µ 2 ) =d(σ, G [0,t 1 ] µ 1 ) + d(σ, F [t 1 ,t 2 ] µ 2 )+ d(σ, G [t 1 ,t min ] (µ 1 ∨ µ 2 )) (6.3)
where t min = min(arg min

t∈[t 1 ,t 2 ]
d e (σ@t, S µ 2 )) is the earliest time corresponding to the closest point between σ and region µ 2 .

These distances are illustrated on Figure 6.2.

d(σ, F [t 1 ,t 2 ] µ) is a piecewise function.
The first condition t last < t 1 and σ@t last ̸ ∈ S µ means that the last transition/jump of the trajectory occurs before t 1 , and the last state of the trajectory does not belong to S µ . In this case, the satisfiability distance computes the Euclidean distance between (t last , σ@t last ) and S µ [t 1 ,t 2 ] . This is related to the choice of ranking trajectories when no reactions occur before t 1 .

Indeed, the efficiency (but not correctness) of our new method for the estimation of Definition 4.3.2 relies on a particular aspect of satisfiability distances, which is not included in the definition. It should compare trajectories fairly by how far they are from satisfying the formula. Suppose σ 1 , σ 2 ∈ P ath(M). If σ 1 is further to fulfil φ than σ 2 , then the distance should verify d(σ 1 , φ) > d(σ 2 , φ) for better efficiency.

An example of a trajectory is depicted in Figure 6.1. d 1 is the satisfiability distance of the violet trajectory because no reactions occur after t 1 . d 2 would be the computed distance by the second sub-function in Equation 6.1. The closer t last is to t 1 , the smaller d 1 . But d 1 is consistently greater than d 2 . This makes a good ranking of trajectories in which no reaction occurs after t, t < t 1 . This is particularly the case when a trajectory reaches an absorbing state before t 1 : it implies no possible evolution of the system after t last that could make the trajectory fulfil the formula.

In this case, the further t last is from t 1 , the further the trajectory is from satisfying the formula.

Satisfiability distances

Thus, by definition, d e (σ@t * , S µ ) = min

s ′ ∈Sµ d i=1 (σ@t * [i] -s ′ [i]) 2 . As s * ∈ S µ , d e (σ@t * , S µ ) = 0. In fine, d(σ, F [t 1 ,t 2 ] µ) = 0. ⇐ ⇐ Suppose that d(σ, F [t 1 ,t 2 ] µ) = 0. Let us prove σ |= F [t 1 ,t 2 ] µ.
As d(σ, F [t 1 ,t 2 ] µ) = 0, there are two possible cases.

First case: t last < t 1 and σ@t last ̸ ∈ S µ . d e ((t last , σ@t last ), S µ [t 1 ,t 2 ] ) = 0 implies

t last = t 1 : impossible. Second case: ∃t * ∈ [t 1 , t 2 ]/d e (σ@t * , S µ ) = 0.
By definition, d e (σ@t * , S µ ) = min

s ′ ∈Sµ d i=1 (σ@t * [i] -s ′ [i]) 2 .
There exists s * ∈ S µ with

d i=1 (σ@t * [i] -s * [i]) 2 = 0, which implies σ@t * = s * . As t * ∈ [t 1 , t 2 ] and σ@t * = s * ∈ S µ , the trajectory σ verifies F [t 1 ,t 2 ] µ. Lemma 6.2.2
For an elementary proposition µ and an MITL temporal formula

G [t 1 ,t 2 ] µ, σ |= φ ⇐⇒ d(σ, φ) = 0 Proof. ⇒ ⇒ Suppose that σ |= G [t 1 ,t 2 ] µ. Let us prove d(σ, G [t 1 ,t 2 ] µ) = 0. As σ |= G [t 1 ,t 2 ] µ, ∀t ∈ [t 1 , t 2 ], σ@t ∈ S µ . We have ∀s ′ ∈ S µ , d e (s ′ , S µ ) = 0. So ∀t ∈ [t 1 , t 2 ], d e (σ@t, S µ ) = 0, hence t 2 t 1 d e (σ@t, S µ )dt = 0. In conclusion, d(σ, G [t 1 ,t 2 ] µ) = 0. ⇐ ⇐ Suppose that d(σ, G [t 1 ,t 2 ] µ) = 0. Let us prove σ |= G [t 1 ,t 2 ] µ.
We have:

t 2
t 1 d e (σ@t, S µ )dt = 0. As t → d e (σ@t, S µ ) is a non-negative continuous function and its integral equals 0, then this function is the null function over

[t 1 , t 2 ]. So, ∀t ∈ [t 1 , t 2 ], d e (σ@t, S µ ) = 0. This means ∀t ∈ [t 1 , t 2 ], σ@t ∈ S µ . In other words, σ |= G [t 1 ,t 2 ] µ.
Now we can tackle the proof of Proposition 6.2.1.

Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs 119 Proof.

d(σ, F [t 1 ,t 2 ] µ i ) = 0 ⇔ min i d(σ, F [t 1 ,t 2 ] µ i ) = 0 ⇔ ∃i/d(σ, F [t 1 ,t 2 ] µ i ) = 0 ⇔ ∃i/σ |= F [t 1 ,t 2 ] µ i ⇔ σ |= F [t 1 ,t 2 ] µ i Result 6.2.1 Def. of disjunction d(σ, G [t 1 ,t 2 ] µ i ) = 0 ⇔ min i d(σ, G [t 1 ,t 2 ] µ i ) = 0 ⇔ ∃i/d(σ, G [t 1 ,t 2 ] µ i ) = 0 ⇔ ∃i/σ |= G [t 1 ,t 2 ] µ i ⇔ σ |= G [t 1 ,t 2 ] µ i Result 6.2.2

Def. of disjunction

For the Until formula, the result follows from

σ |= µ 1 U [t 1 ,t 2 ] µ 2 ⇔ σ |= G [0,t 1 ] µ 1 ∧ G [t 1 ,t min ] µ 1 ∧F [t min ,t 2 ]
µ 2 with the proof of the distances for eventual and global regions.

Linear Hybrid Automata to compute satisfiability distances

Based on the HASL formalism (Section 4.4), we discuss the definitions of hybrid automata to measure the satisfiability distances with synchronised simulations (Section 4.4.2). The automata presented here refer to temporal formulae built on top of a generic mono-dimensional (elementary) proposition µ = x 1 ≤ x O ≤ x 2 (where x O denotes the population of an observable quantity O of a CTMC model M, and

x 1 < x 2 ∈ N).
Each synchronised simulation of a trajectory is associated with the HASL trajectory expression Last(d). The last value of the computed distance variable d is kept.

The Linear Hybrid Automata update the distance on the fly of the simulation. Remark 6.3.1 (Dashed red edge in the automata)

In each automaton, the dashed red edge corresponds to an asynchronous edge that is fired if the current state of the CTMC is absorbing, even if it does not fulfil the condition based on the time t. ---------→, where E(x) is the total exit rate of the state. This implies the definition of a variable test_abs that equal true if no more reaction can occur, which is updated at each 6.3. Linear Hybrid Automata to compute satisfiability distances reaction occurence by the total exit rate. It allows for the acceptance by the automaton of any trajectory, even if it has reached an absorbing state. This variable will not be mentioned in the presentation of the three automata.

More precisely

Distance automaton A F

Automaton A F (Figure 6.3) is designed to measure the satisfiability distance (Equation 6.1) of trajectories of a CTMC model M from the time-bounded region associated with Three variables are defined:

F [t 1 ,t 2 ] (x 1 ≤ x O ≤ x 2 ), i.e. the time-bounded satisfiability region S [t 1 ,t 2 ] (x 1 ≤x O ≤x 2 ) . l 0 ṫ:1 l 1 ṫ:1 l 2 ṫ:1 l 3 ṫ:1 ♯,⊤,{n:=x O ,d:=∞} ♯,t≥t 1 ∧d=0, ∅ ♯,t≥t 2 ∧(n>x 2 ∨x 1 >n), {∅} ♯,t≤t 1 ∧(x 1 >n∨n>x 2 ), {d:=min( √ (t-t 1 ) 2 +(n-x 2 ) 2 , √ (t-t 1 ) 2 +(n-x 1 ) 2 )} ♯,(x 1 ≤n≤x 2 ),{d:=0} ♯,t≥t 1 ∧(x 1 >n∨n>x 2 ), {d:=min(d,min(|n-x 1 |,|n-x 2 |))} ALL,⊤,{n:=x O } t≥t 2 ∨(t≥t 1 ∧d=0),{∅}
• d, the satisfiability distance,

• t, the current time,

• n, the population of the observed species O.

The initial location is l 0 , and the final location is l 2 . First, the distance is initialised to d := ∞ (l 0 → l 1 ), and the initial value of the observed species is stored

in n := x O .
Once in l 1 , there are two ways of reaching the final location l 2 . If the distance d is zero and time exceeds t 1 , the trajectory has entered the time-bounded satisfiability region. Then, the transition l 1 t≤t 1 ∧d=0,∅ -------→ l 2 fires. Otherwise, if the time exceeds t 2 , the second transition to l 2 is traversed without updating the variables.

If, while in l 1 , the trajectory has not entered the time-bounded satisfiability region yet, the distance d must be updated. It deals with three autonomous transitions l 1 → l 3 .

On the one hand, in case x O has entered [x 1 , x 2 ], the distance is set to d := 0 (l 1 (x 1 ≤n≤x 2 ),{d:=0} ------------→ l 3 ). Indeed, since CTMC trajectories are stepwise functions, if the next reaction occurs at time t > t 1 , the current trajectory has at least one point within the time-bounded satisfiability region.

On the other hand, if ------------------------------→ l 3 fires: it computes the distance between the trajectory and the nearest corner of the eventual region. It corresponds to the first sub-function of the piecewise Equation 6.1.

x O has not entered [x 1 , x 2 ], d is computed with two different equations. First, if t ≤ t 1 , l 1 ♯,(x 1 >n∨n>x 2 )∧(t≤t 1 ), {d:=min( √ (t-t 1 ) 2 +(n-x 2 ) 2 , √ (t-t 1 ) 2 +(n-x 1 ) 2 )} --
Second, if t ≥ t 1 , the third transition fires. d is computed as the minimum distance between d and the distance of n to the bounds x 1 and x 2 , which corresponds to the second sub-function of the piecewise Equation 6.1.

In l 3 , two edges can be traversed. First, the asynchronous transition brings the automaton to the final location l 2 if either the distance is zero within [t 1 , t 2 ], time exceeds t 2 or the trajectory reaches an absorbing state. In addition, the transition l 3 ALL,true,{n:=x O } -----------→ l 1 is traversed whenever a new reaction occurs, which brings the automaton back to l 1 so that the distance can be updated accordingly.

Remark 6.3.2 (Left-closed constraints on transition predicates)

In the definition of an LHA, transition predicates of autonomous transitions have to be left closed linear constraints (Definition 4.4.2). But, for example t ≥ t 2 ∧ (n > x 2 ∨ x 1 > n) is not a left-closed constraint. Thus, the autonomous transition l 1 ♯,t≥t 1 ∧(x 1 >n∨n>x 2 ),d=...

-----------------→ l 3 defines, in fact, two autonomous transitions:

l 1 ♯,t≥t 1 ∧x 1 +1≥n,d=...
-------------→ l 3 and l 1 ♯,t≥t 1 ∧n≥x 2 +1,d=...

-------------→ l 3 as n, x 1 and x 2 are integers.

Script 7 (A F LHA in Cosmos)

An example of automaton A F in Cosmos is available at code/chap5/cosmos/ER/ dist_F.lha. Automaton A G (Figure 6.4) is designed to measure the satisfiability distance (Equation 6.2) of trajectories of a CTMC model M from the time-bounded region associated

Distance automaton A

G l 0 ṫ:1 l 1 ṫ:1 l 3 ṫ:1 ṫ′ :0 l 4 ṫ:1 ṫ′ :1 l 2 ṫ:1 ♯,tru e,{n :=x O , d:= 0,in :=tr ue} ♯,in∧t 1 ≤t≤t 2 ∧(x 1 ≤n≤x 2 ), {t ′ :=0} ♯,t≤t 1 ∧(x 1 ≤n≤x 2 ) {d:=0,in:=false} ♯,¬in∧t 1 ≤t≤t 2 ∧(x 1 ≤n≤x 2 ) {d:=d•(t-t 1 ),t ′ :=0} ♯,t≤t 1 ∧(n<x 1 ∨n>x 2 ), {d:=min(|x 1 -n|,|x 2 -n|),in:=false} ♯ ,¬ in ∧ t 1 ≤ t ≤ t 2 ∧ ( n < x 1 ∨ n > x 2 ) , { d + = d • ( t -t 1 ) } ♯ ,i n ∧ t1 ≤ t ≤ t2 ∧ ( n < x1 ∨ n > x2 ) , ∅ ALL,true,{n:=x O } ♯ , i n ∧ t ≥ t 2 { d = d • ( t 2 - t 1 ) } ♯ , ¬ i n ∧ t ≥ t 2 , { ∅ } A L L ,t ru e, { d + = t ′ •m in (| x1 -n |, |x2 -n |) )} , { t ′ := 0 ,n := x O ,i n := tr ue } ♯,t≥t 2 ,{d+=t ′ •min(|x 1 -n|,|x 2 -n|))}
with

G [t 1 ,t 2 ] (x 1 ≤ x O ≤ x 2 ), i.e.
(t ′′ -t 1 ) * min(|n -x 1 |, |n -x 2 |) instead of (t ′′ -t ′ ) * min(|n -x 1 |, |n -x 2 |).
After the initialisation of the variables (l 0 → l 1 ), analysis begins in l 1 . For jumps occurring before t 1 , we distinguish two cases. If σ@t ∈ [x 1 , x 2 ], the distance is set to zero (l 1 → l 3 top arc). Otherwise, d is the distance of σ@t from [x 1 , x 2 ] (l 1 → l 3 midway arc). Indeed, if, for example, the next jump of σ happens at t > t 2 , then the final distance is given by d • (t 2 -t 1 ) (l 1 → l 2 bottom arc).

For jumps occurring within t

∈ [t 1 , t 2 ], if σ@t ̸ ∈ [x 1 , x 2 ] (sequence l 1 → l 4 → l 1 ),
the distance is incremented by the surface contoured with the trajectory segment (of duration t ′ ) laying outside [x 1 , x 2 ], and the closest border of

[x 1 , x 2 ]. The distance is left unchanged if σ@t ∈ [x 1 , x 2 ] (sequence l 1 → l 3 → l 1 ). Script 8 (A G LHA in Cosmos)
An example of automaton A G in Cosmos is available at code/chap5/cosmos/ER/ Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs 123 dist_G.lha.

Distance automaton A G∧F

Automaton A G∧F (Figure 6.5) is designed to compute the distance of trajectories from a sequence of time-bounded regions consisting of a formula

G [t 1 ,t 2 ] (x 1 ≤ x O ≤ x 2 )
(related to an observed quantity O) followed by a formula

F [t 3 ,t 4 ] (x 3 ≤ x O ′ ≤ x 4 )
(related to an observed quantity O ′ ). This automaton is associated with the MITL

formula φ = G [t 1 ,t 2 ] (x 1 ≤ x O ≤ x 2 ) ∧ F [t 3 ,t 4 ] (x 3 ≤ x O ′ ≤ x 4 ). We assume t 2 ≤ t 3 , i.e.
the global region precedes the eventual region, while x 1 , x 2 , x 3 , x 4 ∈ N, and x O , resp.

x O ′ , denotes the population of species O, resp. O ′ .

In particular, it deals with the satisfiabily distance of conditional time-bounded

reachability because G [0,s] (x 1 ≤ x O ≤ x 2 ) ∧ F [s,t] (x 3 ≤ x O ′ ≤ x 4 ) = (x 1 ≤ x O ≤ x 2 )U [s,t] (x 3 ≤ x O ′ ≤ x 4 ) when µ = (x 1 ≤ x O ≤ x 2 )∨(x 3 ≤ x O ′ ≤ x 4
) is verified for all states of any possible trajectory of the CTMC (see Remark 6.2.1).

We combine the automata A G and A F by linking them with an autonomous transition l 2G → l 1F . d stores the distance for the G region, while d ′ stores the distance for the F region. The sum of the two distances is stored in d at the end of the simulation.

6.4. Automaton-ABC algorithm with LHA satisfiability distances Script 9 (A G∧F LHA in Cosmos)

l 0G ṫ:1 l 1G ṫ:1 l 3G ṫ:1 ṫ′ :0 l 4G ṫ:1 ṫ′ :1 l 2G ṫ:1 l 1F ṫ:1 l 2F ṫ:1 l 3F ṫ:1 ♯,tru e,{n :=x O , d:= 0,in :=tr ue} ♯,in∧t 1 ≤t≤t 2 ∧(x 1 ≤n≤x 2 ), {t ′ :=0} ♯,t≤t 1 ∧(x 1 ≤n≤x 2 ) {d:=0,in:= false} ♯,¬in∧t 1 ≤t≤t 2 ∧(x 1 ≤n≤x 2 ) {d:=d•(t-t 1 ),t ′ :=0} ♯,t≤t 1 ∧(n<x 1 ∨n>x 2 ), {d:=min(|x 1 -n|,|x 2 -n|),in:=false} ♯ , ¬ in ∧ t 1 ≤ t ≤ t 2 ∧ ( n < x 1 ∨ n > x 2 ) , { d + = d • ( t - t 1 ) } ♯ ,i n ∧ t 1 ≤ t ≤ t 2 ∧ ( n < x 1 ∨ n > x 2 ) , ∅ i n ∧ t ≥ t 2 { d = d • ( t 2 - t 1 ) } ¬ i n ∧ t ≥ t 2 , { ∅ } ALL,true,{n:=x O } A L L ,t ru e ,{ d ′ + = t ′ • m in ( | x 1 -n | , | x 2 -n | ) ) } , { t ′ := 0 , n := x O , in := tr u e } ♯,t≥t 2 ,{d+=t ′ •min(|x 1 -n|,|x 2 -n|))} ♯, tru e,{ n: = x O ′ ,d ′ := ∞ } ♯,t≥t 1 ∧d ′ =0, {∅} ♯,t≥t 2 ∧(n>x 2 ∨x 1 >n), {d=d+d ′ } ♯, t≤t 1 ∧(x 1 >n∨n>x 2 ), {d ′ :=min( √ (t-t 1 ) 2 +(n-x 2 ) 2 , √ (t-t 1 ) 2 +(n-x 1 ) 2 )} ♯,(x 1 ≤n≤x 2 ),{d ′ :=0} ♯,t≥t 1 ∧(x 1 >n∨n>x 2 ), {d ′ :=min(d ′ ,min(|n-x 1 |,|n-x 2 |))} ALL,true,{n:=x O ′ } t≥t 2 ∨(t≥t 1 ∧d=0),{d=d+ d ′ }
An example of automaton A G∧F in Cosmos is available at code/chap5/cosmos/ER/ dist_G_F.lha.

Automaton-ABC algorithm with LHA satisfiability distances

In this part, we define the automaton-ABC framework with satisfiability distances to estimate the satisfaction probability function (Section 6.1). The framework aims to first find the subset of parameters of a parametric CTMC model (M θ ) θ∈Θ for which a reachability formula φ can be satisfied by efficiently exploring the whole parameter space, and then estimate the satisfaction function probability. The intuition behind this method is that the exploration of the parameter space can be driven efficiently by taking into account the notion of distance of a trajectory σ ∈ P ath(M θ ) from the timebounded satisfiability region corresponding to φ (see Section 6.2). In this formulation of the ABC algorithm, the estimation of the posterior distribution (π φ-ABC ) is no longer computed as a limit approximation (i.e. lim ϵ→0 π ϵ ABC (.|y exp )), as with classical ABC, but rather as an estimation of the exact ABC distribution, since trajectories are accepted exclusively if their distance to the satisfiability region is zero.

Simple ABC with satisfiability distance.

We consider a parametric CTMC (M θ ) θ∈Θ , and we define a prior distribution π over Θ. In Algorithm 17, we propose a modified version of the general automaton-ABC Algorithm 15 adapted to satisfiability distances. At each iteration, we draw a parameter θ ′ from the prior π(.), we simulate a path σ ′ from the CTMC M θ ′ synchronised with an LHA A φ (Section 4.4.2) and accept θ ′ if the distance from φ is d(σ ′ , φ) = 0. (i.e. if σ ′ |= φ by Proposition 6.2.1).

Algorithm 17 Automaton-ABC with A φ automaton Require: (M θ ) θ∈Θ a pCTMC, π(.) prior, A φ automaton distance for MITL formula φ, N : number of particles Ensure: (θ (i) ) 1≤i≤N drawn from π φ-ABC

for i = 1 : N do repeat θ ′ ∼ π(.) d(σ ′ , φ) ∼ (Last(d), A φ ) × M θ ′ until d(σ ′ , φ) = 0 θ (i) ← θ ′ end for
The Theorem 6.4.1 below links the Algorithm 17 with the satisfaction probability function.

Theorem 6.4.1

Let (M θ ) θ∈Θ be a parametric CTMC, φ an MITL formula, and the associated satisfaction probability function θ → P r(φ; M θ ). Given a prior distribution π over the parameter set Θ, the (θ (i) ) 1≤i≤N sampled by the Algorithm 17 are drawn from a density π φ-ABC :

π φ-ABC (θ) = P r(φ; M θ ) π(θ) C
where C ∈ R ≥0 is a positive constant.

6.4. Automaton-ABC algorithm with LHA satisfiability distances

Proof. Let C φ M θ = {σ ∈ P ath(M θ )/σ |= φ}. Then P r M θ (C φ M θ ) = P r(φ; M θ ).
ABC is a reformulation of the accept-reject algorithm (see Section 3.2.2 and Algorithm 11). Thus, the samples (θ (i) , σ (i) ) 1≤i≤N from Algorithm 17 are drawn from a density π φ-ABC :

π φ-ABC (θ (i) , σ (i) ) ∝ 1 d(.,φ)=0 (σ (i) ) σ (i) |=φ p M θ (i) (σ (i) ) σ (i) ∼M θ (i) π(θ (i) )
prior where p M θ (i) is the density related to the CTMC M θ (i) with respect to a measure µ.

As i) ). One can obtain the marginal distribution of θ by integration over the whole set of trajectories P ath(M θ ):

σ (i) ∈ C φ M θ (i) ⇔ d(σ (i) , φ) = 0, 1 d(.,φ)=0 (σ (i) ) = 1 C ϕ M θ (i) (σ (
π φ-ABC (θ) ∝ σ∈P ath(M θ ) 1 C ϕ M θ (σ)p M θ (σ)π(θ)dµ ∝ π(θ) σ∈C ϕ M θ p M θ (σ)dµ ∝ P r M θ (C φ M θ )π(θ) density of M θ
As P r M θ (C φ M θ ) = P r(φ; M θ ), we can conlude that there exists C ∈ R ≥0 so that π φ-ABC (θ) = P r(φ; M θ ) π (θ) C .

This result transforms the regression of the smooth function f φ into the density estimation of π φ-ABC . First, each parameter sampled from the π φ-ABC gives information because it produces a simulation that verifies φ. Also, as π φ-ABC is a probability density function, the relative position of each parameter to the other sampled parameters gives much information about the satisfaction probability function. The denser in sampled parameters a subset of parameters space, the higher the satisfaction probability function over the subset.

Following the same approach, we formulate a comparable version for the ABC-SMC Algorithm 18, which leads to a smaller runtime than Algorithm 17, keeping the guarantee of Theorem 6.4.1. Indeed, Algorithm 17 do not really exploit the continuity of the distance. We only accept/reject trajectories depending on whether their distance is zero, i.e. if they satisfy φ. On the contrary, with Algorithm 18, we use the satisfiability distance to rank paths and accept the parameters whose corresponding paths are closer to the satisfiability regions (better ranked) than others, even if they do not necessarily satisfy φ. This new feature can lead to faster convergence of the algorithm corresponding to a faster exploration of the parameter space.

Algorithm 18 has the same inputs as Algorithm 17, plus a kernel distribution K and a hyperparameter α ∈]0, 1[ representing how fast the tolerance ϵ decreases with the iterations.

It works as follows. Initially, N parameters/particles (θ

Estimation of the satisfaction probability function

Based on the samples (θ (i) ) 1≤i≤N (Algorithm 17 or Algorithm 18), we can estimate the satisfaction probability function thanks to Theorem 6.4.1. This procedure is twofold: estimation of the ABC posterior density and estimation of the constant C.

Estimation of the π φ-ABC posterior distribution

In the applications, we estimate the ABC posterior π φ-ABC based on the samples (θ (i) ) 1≤i≤N with kernel density estimation (Section 3.4

): ∀θ ∈ R, πh (θ) = 1 N d i=1 K h (θ, θ (i) ).
Two kernels are used: Gaussian and Beta kernels. The last one is helpful when we have to estimate densities over bounded supports with positive probabilities on the boundaries. The optimal bandwidth is obtained by Least Squares Cross-Validation.

6.4. Automaton-ABC algorithm with LHA satisfiability distances Algorithm 18 Automaton-ABC Sequential Monte Carlo with A φ automaton Require: (M θ ) θ∈Θ a pCTMC, π prior, A φ distance automaton for MITL formula φ, N : number of particles, α ∈ (0, 1), K kernel distribution Ensure: (ω (i) , θ (i) ) 1≤i≤N weighted samples drawn from π φ-ABC

θ (i) 0 ∼ π, i ∈ 1, . . . , N d i ∼ (Last(d), A φ ) × M θ (i) 0 , i ∈ 1, . . . , N ϵ ← quantile(α, (d i ) 1≤i≤N ) (ω (i) 0 ) 1≤i≤N ← 1 N m ← 1 while ϵ > 0 do for i = 1 : N do repeat Take θ ′ from (θ (j)
m-1 ) 1≤j≤N with probabilities (ω

(j) m-1 ) 1≤j≤N θ (i) m ∼ K(.|θ ′ ) d i ∼ (Last(d), A φ ) × M θ (i) m until d i ≤ ϵ ω (i) m ← π ( θ (i) m ) N Σ i ′ =1 ω (i ′ ) m-1 K(θ (i) m |θ (i ′ ) m-1 )
end for Normalise (ω

(i) m ) 1≤i≤N ϵ ← quantile(α, (d i ) 1≤i≤N ) m ← m + 1 end while return (ω (i) m , θ (i) m ) 1≤i≤N 6.5. Applications The set C φ M λ is the set of trajectories of M λ that satisfies φ = F [t 1 ,t 2 ] (x 1 ≤ s ≤ x 2 )
for the CTMC M λ . There are powerful known results on Poisson processes such as:

i. The trajectories are increasing step functions with jumps of value one (i.e has the form 1

t 1 -→ 2 t 2 -→ 3 → . . .). ii. ∀t, v ∈ R ≥0 with v ≤ t, S λ t -S λ,v ∼ P(λ(t -v)), i.e. a Poisson distribution of parameter λ(t -v).
With the property i), having one point in the time-bounded satisfiability region of F [t 1 ,t 2 ] (x 1 ≤ s ≤ x 2 ) is equivalent to having σ@t 1 ≤ x 2 and σ@t 2 ≥ x 1 . Indeed, if σ@t 1 > x 2 , the trajectory is on the top of the region, and if σ@t 2 < x 1 , the trajectory is below the region. Else, the trajectory must traverse [x 1 , x 2 ].

With this remark, P

S λ (C φ M λ ) = P(S λ t 1 ≤ x 2 , S λ t 2 ≥ x 1 ). P(S λ t 1 ≤ x 2 , S λ t 2 ≥ x 1 ) = P(S λ t 2 ≥ x 1 |S λ t 1 ≤ x 2 )P(S λ t 1 ≤ x 2 ) = x 2 k=0 P(S λ t 2 ≥ x 1 |S λ t 1 = k)P(S λ t 1 = k|S λ t 1 ≤ x 2 ) = x 2 k=0 P(S λ t 2 -S λ t 1 ≥ x 1 -k)P(S λ t 1 = k) = x 2 k=0 (1 -P(S λ t 2 -S λ t 1 < x 1 -k))P(S λ t 1 = k) = x 2 k=0 (1 -P(S λ t 2 -S λ t 1 ≤ x 1 -k -1))P(S λ t 1 = k) {S λ t1 = k} k=[[0,x2]] is a partition of {S λ t1 ≤ x 2 }
As S λ t 2 -S λ t 1 ∼ P(λ(t 2 -t 1 ))) and S λ t 1 ∼ P(λt 1 ), P S λ (C φ M λ ) can be expressed analytically. As P S λ = P r M λ (Proposition 2.2.1), the satisfaction probability function

λ → P r M λ (σ |= F [t 1 ,t 2 ] (x 1 ≤ s ≤ x 2 )
) is expressed analytically. We have solved the time-bounded reachability problem analytically, which offers the possibility to test our method with exact continuous solutions.

The Figure 6.6 below shows the run of the automaton-ABC algorithm with a prior π(.) ∼ U (0, 30), N = 500 particles with x 1 , x 2 = 5.0, 7.0 and t 1 , t 2 = 0.75, 1.0.

The estimated satisfaction probability function with Gaussian kernel estimation of the ABC posterior is illustrated in red, whereas the true satisfaction probability function is in blue. One can see that our estimation suits well the exact satisfaction probability function. In blue: the analytical satisfaction probability function.

Enzymatic reaction system

Model

We consider the Enzymatic Reaction system (Michaelis-Mentens kinetics [START_REF] Michaelis | The Original Michaelis Constant: Translation of the 1913 Michaelis-Menten Paper[END_REF]) described in Equations 6.4 in which a substrate species S is converted into a product P through the mediation of an enzyme E. The dynamics depend on the reaction rate constants that define a parameter vector θ = {k 1 , k 2 , k 3 }. We thus consider the underlying parametric CTMC (M) θ∈[0,100] 3 . The initial state is (E 0 , S 0 , ES 0 , P 0 ) = (100, 100, 0, 0). The dynamic of the ER system (Figure 6.7) is such that the totality of the substrate (initially S 0 = 100) is converted into the product at a speed dependent on parameters θ. With θ = (1, 1, 1), the totality of S is converted before time = 5, whereas with a tenfold speed reduction in the formation of the ES complex and 6.5. Applications synthesis of P (i.e. θ = (0.1, 1, 0.1)), we have that only about 30% of S has been converted at time = 5.

R 1 : E + S k 1 -→ ES R 2 : ES k 2 -→ E + S R 3 : ES k 3 -→ E + P (6.4)

Experimental setting

In the following, we run a few experiments on the ER system. We consider six time-bounded reachability formulae and their associated time-bounded satisfiability regions:

• φ 1 : F [0.025,0.05] (50 ≤ P ≤ 75) associated with region T R1,

• φ 2 : F [0.05,0.075] (50 ≤ P ≤ 75) associated with region T R2,

• φ 3 : F [0.05,0.075] (25 ≤ P ≤ 50) associated with region T R3, Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs

133

• φ 4 : F [8.0,10.0] (5 ≤ P ≤ 15) associated with region T R4,

• φ 5 : G [0.0,0.8] (50 ≤ E ≤ 100) associated with region T R5,

• φ 6 : G [0.0,0.8] (50 ≤ E ≤ 100) ∧ F [0.8,0.9] (30 ≤ P ≤ 100) associated with region T R6.

φ 1 , φ 2 , φ 3 correspond to formulae that are probably fulfilled when k 3 varies over [0, 100] (with k 1 = k 2 = 1), whereas φ 4 , φ 5 , φ 6 are formulae for which only a small subset of k 1 , k 2 ∈ [0, 100] 2 (with k 3 = 1) can produce trajectories that fulfil these formulae. Such intuition is confirmed by plots in Figure 6.9, which depicts the average value of the distance of trajectories from time regions T R1, T R2 and T R3 as a function of k 3 . We use the Statistical Model Checker Cosmos [START_REF] Ballarini | Applications of an expressive statistical model checking approach to the analysis of genetic circuits[END_REF] with specific instances of A F i.e. A φ 1 , A φ 2 and A φ 3 . We observe that for example the measured distance from region T R1 monotonically decreases as k 3 increases and cancels for k 3 ≥ 30, while the distance from region T R3 is zero when 10 ≤ k 3 ≤ 15 whereas it grows as we increase k 3 .

Test of LHA distances

Satisfaction probability function estimation

We apply the automaton-ABC Algorithm 18 to the ER parametric CTMC defined over [0, 100] 3 , with formulae described in Section 6.5.2.

T R1, T R2 and T R3 correspond to the one-dimensional experiments: only k 3 varies over [0, 100] (k 1 = k 2 = 1.0 are fixed), a uniform prior π(.) ∼ U (0, 100) is set.

T R4, T R5 and T R6 correspond to the two-dimensional experiments: k 1 and k 2 varies over [0, 100] (k 3 = 1.0 is fixed), a uniform prior π(.) ∼ U (0, 100) over each parameter is set.

T R1 T R2 T R3

F [0.025,0.05] (50 ≤ P ≤ 75) F [0.05,0.075] (50 ≤ P ≤ 75) F [0.05,0.075] (25 ≤ P ≤ 50) 6.5. Applications entering T R4, or never leaving T R5, which means that the algorithm managed to catch the correlation between the parameters. This is intuitively correct in both cases. In fact, T R4 corresponds to a very low synthesis of P , which is not compatible with fast creation of the ES complex (i.e. only very small k 1 are not ruled out), and even the compensation effect obtained by fast decomplexation (i.e. large k 2 ) will not suffice for the different trajectories to stay in T R4.

One can notice that the estimated satisfaction probability function of region T R4

(most-left bottom picture of Figure 6.11) has much probability mass around (0, 0).

At first glance, one can conclude in a problem of bias of the kernel density estimator since the satisfaction probability function estimated by Statistical Model Checking (second row, first column picture) does not show the same shape in this area. But, if we execute Statistical MC with a refined grid around zero (Figure 6.12), one can see surprisingly that probability values are high around (0, 0). This behaviour was not expected and automaton-ABC algorithm permits to discover this small area of high satisfaction probability that was not caught by Statistical Model Checking over the grid of 20 points per axis.

Similarly to experiment T R4, T R5 limits the speed of the initial decrease of E (which is initially E 0 = 100), to 50 within t ≤ 0.8. Again, this behaviour is only compatible with a slow ES complexation and cannot be compensated by fast decomplexation. The T R6 experiment caught an even more important correlation between parameters. In fact, the support of T R6 posterior is contained within the support of T R5 posterior. Indeed, a trajectory that verifies φ 6 also verifies φ 5 .

Applications

Intracellular viral infection

We consider a model of cell viral infection [START_REF] Haseltine | Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics[END_REF] described by Equations 6.5.4.

N + T k 1 X T cn -----→ G + T T k 4 X T ---→∅ N + G k 2 X G cn -----→ T S k 5 X S ---→∅ N + A + T k 3 X T cnca ------→ S + T G + S k 6 X G X S -----→V
N represents the nucleotides and A the amino acids. G is the genomic nucleic acids, T the template nucleic acids, S the viral structural protein and V the secreted virus. In this model, we assume nucleotides and amino acids have constant concentrations c n and c a : we suppose these species are in large number.

Satisfaction probability function estimation A 2D experiment is considered

with the following logical property: φ = G [0,50] (G ≤ 10) ∧ F [50,200] (G > 100). This property expresses that the species G remains stable and low within time [0, 50], and then a burst of speed in the creation of G occurs within [50, 100], which gives the necessary material for the creation of the virus V .

We vary the concentrations of nucleotides and amino acids: c n ∼ U (0.6, 1.1)

and c a ∼ U (0.5, 2.0). Figure 6.14 reports the results of the experiment. Considering the end time of the formula, each simulation of this model is computationally expensive. But our methodology still allows a proper run of the experiment, and we get a reasonable approximation of the satisfaction function (it required about 7 • 10 3 simulations).

About the implementation of the automaton-ABC method

The implementation of automaton-ABC methods led to the development of a Julia package. It includes both simulation and synchronised simulation of CTMC, easy modeling of CTMCs based on Chemical Reaction Network, and ABC related methods.

Many details are provided in Appendix A.

ABC related algorithms are distributed and the above applications were performed using HPC resources from the "Mésocentre" computing center of Centrale-Supélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région

Île-de-France (http://mesocentre.centralesupelec.fr/). Tables 6.1 and 6.2 show the performance results for the whole set of experiments.

We omit the Kernel Density Estimation. It only depends on the number of particles N , the length of the grid, and the choice of the kernel we set. Still, the computational time of Least Squares Cross-Validation can be prohibitive when the kernel is a multivariate beta kernel with a high number of particles (N ≥ 1000) and grid search (1000 is already high even if the computation of LSCV is distributed).

When the number of jobs is 120, the run was distributed on the Mesocentre HPC cluster. Otherwise, it was run on a Dell XPS 9370 with CPU Intel i7 8550U @ 1.8GHz

x 8 cores.

Our tool for automaton-ABC in Julia is very efficient because, in our applications, any experiment run is rarely higher than one minute. For example, the experiment about the region TR4 of the ER model performs 256000 simulations within a minute over the cluster (without counting the other computations of the ABC algorithm), knowing that each simulation that reaches T R4 is about 2000 steps of Stochastic Simulation Algorithm 1. In Appendix A, we compare the performance of our tool with state-of-the-art softwares.

The SIR 2D experiment has a higher computational time than the 1D experiment, even if it is run with 120 jobs. This is explained by the fact that creating and dealing with many jobs has a higher computational cost when the execution time per job is low. It is the case here: less than 7 seconds of computation has to be distributed over 120 jobs, which is not helpful. The viral infection model has a higher computational cost per simulation because it is a much more complex model than the two others. In return, the ABC method does not have the same statistical guarantees as Smoothed MC. Indeed, we do not assume any specific form for the density π φ-ABC :

we minimise Least Squares Cross-Validation criterion to have the better trade-off between bias and variance over the N = 1000 particles in kernel density estimation.

To run the 2D experiments, we have adapted the available Smoothed MC code, but preliminary tests seem to indicate a prohibitive computational time. Indeed, convergence for 2D experiments such as T R4 required many more trajectories simulation (∼ 2 • 10 5 ). The end time of simulations is much higher than the 1D experiments (t 2 = 10.0 ≫ 0.075).

The 2D experiments on the ER system show that our method gives an efficient way to identify the region of the parameter space where the satisfaction function is positive. After such an exploration is performed, one can use either kernel density estimation procedure, or use regression function methods such as Smoothed MC, that trade higher computational cost for better statistical guarantees.

Note that since the main computational cost of the ABC approach is the simulation of trajectories, our method is well suited for distributed computing, which is not the case for all regression methods (e.g. Smoothed MC is, by nature, sequential).

Discussion

6.7.1 About the distance of automaton A F before t 1

In Section 6.2, we detailed a particular aspect of the automaton: if the entire system does not evolve between a time t ′ < t 1 until t 2 , the distance computed by A F equal the first part of Equation 6.1. However, with this behaviour, if the system evolves but the number of species x O remains the same, the distance will be updated. For example, in Figure 6.1, if reactions occur between [t 1 , t 2 ] = [8.0, 10.0] that does not affect the copy number of P , the output of A F is d 2 (the smaller one). It can be problematical if there are reactions completely uncorrelated to the others that can still occur (for example, adding a reaction ∅ → A to the ER system).

If one wants to change this behaviour, a possible solution is to add an extra variable n ′ that stores the second to last value of the number of the studied species O. Then, the edges l 1 → l 3 must contain a new predicate n ̸ = n ′ to know whether the automaton should update the distance or not. This approach increases the computational time of the synchronised simulation. This thesis addressed several tasks of statistical inference and verification for Continuous-Time Markov Chains described by Chemical Reaction Networks. We have endeavoured to prove the efficiency of the synchronised simulation of a CTMC with a Linear Hybrid Automaton within the ABC algorithm. This new method called automaton-ABC has been applied to several problems from Systems Biology.

First, we gave an overview of the literature for Bayesian statistical inference and statistical model checking for parametric CTMCs. We introduced the notations from both fields, detailed what makes the considered problems complex, and discussed state-of-the-art techniques.

Second, we used automaton-ABC for statistical inference of oscillatory genetic networks. The nature of observations concerns the oscillations. They are described by sentences such as "we observe N oscillations between values L and H with mean period duration tp (obs) ". This qualitative behaviour is measured by an LHA A per [START_REF] Ballarini | Applications of an expressive statistical model checking approach to the analysis of genetic circuits[END_REF], and we define a distance over such property of an oscillatory trajectory. The use of automaton-ABC over two oscillatory models (doping threeway oscillator and repressilator) allowed to show the parameter space subsets where specific kind of oscillations are the most probable.

Third, we used automaton-ABC to gain computational time in classical inference. HASL formalism offers a framework that tailors the simulation according to a Finally, since these methods mix several concepts from both fields, we have tried to make their use accessible and efficient. We chose to use Julia because of its high-level syntax, efficiency, and great possibility of metaprogramming and code generation. It led to a Julia package that is much detailed in Appendix A. All the Julia scripts related to the results are documented within the Chapters, and the Scripts List above the introduction.

Limits and Perspectives

Scope of our work

This work focused on CTMCs defined by Chemical Reaction Networks. Indeed, they quickly challenge statistical inference or verification problems. However, the methods are easily generalisable to generalised semi-Markov processes (which corresponds to Generalised Stochastic Petri Nets).

Automaton-ABC for statistical inference

In Section 5.3, we gave a new approach to study oscillations in genetic networks whose stochastic dynamics are described by CTMCs. Another qualitative behaviour much observed in biological systems is bistability [START_REF] Tyson | Biological switches and clocks[END_REF]. Bistable models are characterised by the existence of two steady states (and possibly unstable states). However, such stochastic models can produce trajectories that switch between two steady states. Thus, one can imagine an LHA that would characterises this behaviour, for example, by quantifying the time spent in a detected steady state or the duration time of a transition between two steady states, and use this LHA within the automaton-ABC procedure to identify the regions of the parameter space that reproduces this behaviour.

Automaton-ABC for time-bounded reachability

As already mentioned in Section 6.7.2, the method only covers a fragment of MITL.

Further work has to be done to cover all MITL formulae. Instead of trying to design an automaton for each kind of formula, defining operations on HASL trajectory formulae could improve the coverage of MITL.

To illustrate this, we consider the computation of satisfiability distances of MITL formulae φ 1 ∧ φ 2 or φ 1 ∨ φ 2 . We suppose φ 1 and φ 2 are formulae whose satisfiability distances are computed by already defined LHA A φ 1 and A φ 2 with HASL trajectory formulae (Last(d φ 1 ), A φ 1 ) and (Last(d φ 2 ), A φ 2 ). The idea is to define an LHA A φ 1 × A φ 2 that integrates the behaviours of A φ 1 and A φ 2 into one LHA. One can construct this automaton by defining a product over two LHA, whose location set is the cartesian product of A φ 1 and A φ 2 ; variables and flows are the concatenation of those from the two automata. For example, let us consider these two simple LHA:

l 1 l 2 E, γ, U l ′ 1 l ′ 2 E ′ , γ ′ , U ′
Its product could be defined as:

7.1. Limits and Perspectives

(l 1 , l ′ 1 ) begin (l 1 , l ′ 2 ) (l 2 , l ′ 1 ) (l 2 , l ′ 2 ) E, γ, U E, γ, U E ′ , γ ′ , U ′ E ′ , γ ′ , U ′ E ∩ E ′ , γ ∧ γ ′ , U ∪ U ′
With such operations, the distances from both formulae φ 1 and φ 2 would be computed even if simultaneous updates of the two automata variables occur. Indeed, (l 1 , l ′ 1 ) → (l 2 , l ′ 2 ) is fired if the occurred event is in E ∩ E ′ and the predicate γ ∧ γ ′ is fulfilled. Thus, the satisfiability distance of φ 1 ∧ φ 2 could be (Last(d φ 1 + d φ 2 ), A φ 1 × A φ 2 ). For φ 1 ∨ φ 2 , the satisfiability distance could be (Last(d

φ 1 * d φ 2 ), A φ 1 × A φ 2 ).
As [START_REF] Bortolussi | Bayesian statistical parameter synthesis for linear temporal properties of stochastic models[END_REF] did with [START_REF] Bortolussi | Smoothed model checking for uncertain Continuous-Time Markov Chains[END_REF], a statistical parameter synthesis method may be constructed on top of automaton-ABC for time-bounded reachability. The stake is to statistically control the approximation of kernel density estimator based on automaton-ABC samples. For example, one can relax the tolerance by setting ϵ > 0 to have better confidence for parameters that do not verify a formula φ.

In complex CTMC models, exact simulation has a high computational cost. One can use approximate simulation methods such as Tauleaping (Algorithm 3) to reduce the computational cost. The smaller the time resolution τ , the better the simulation approximation. The Multi-Level ABC method [START_REF] Lester | Multi-level Approximate Bayesian Computation[END_REF], evoked in Section 5.1.2, could be a solution to improve the performance of our automaton-ABC procedure.

Implementation

In our implementation detailed in Appendix A, defining a CTMC by a CRN and apply ABC-SMC and automaton-ABC-SMC is easy with the LHA presented in this thesis. However, further work will be performed to improve our tool:

• It is possible but challenging (without help) to design its own LHA. Thus, work is in progress to render the creation of LHA as simple as CTMCs in the package. i. load_model("ER") loads the ER model.

ii. load_automaton("automaton_F") loads a function create_automaton_F that creates an automaton A F described in Figure 6.3 of Chapter 6.

iii. A_F = create_automaton_F(ER, 50.0, 75.0, 0.025, 0.05, :P) creates an automaton A F with constants x 1 = 50, x 2 = 100, t 1 = 0.025, t 2 = 0.05.

iv. σ = simulate(sync_ER) simulates the sync_ER model and stores the simulation in the variable σ.

v. plot(σ) plots the trajectory σ.

Script 13 (Synchronised simulation of the ER model with MarkovProcesses.jl) Here is a detailed run with ϵ = 0.2 and 1000 particies. Abstract: Chemical Reaction Networks (CRN) constitute a formalism used to model biological processes. Under certain assumptions, a Continuous-Time Markov Chain describes its stochastic dynamics: the future state of the system only depends on the current state.

Statistical inference of such CTMCs is complex: likelihood computations are generally intractable. Approximate Bayesian Computation is a recent class of likelihood-free methods for Bayesian inference that allows approximating the posterior distribution with Monte Carlo simulations. It has proven its efficiency in the case of CTMCs. There is a growing interest in the verification (model-checking) of models from Systems Biology to understand the complex molecular interactions within a biological system. Statistical model checking methods have been developed to apply these methods on more complex models. They simulate the model and compute the ratio of simulations that fulfils a property. Recently, Hybrid Automata Stochastic Logic (HASL) has been introduced for the statistical verification of stochastic models. This temporal logic inherently adopts the statistical point of view of model-checking. In this thesis, we focus on statistical inference and verification of CTMCs defined by CRNs.

Our main contribution consists in the new formulation of an Approximate Bayesian Computation procedure combined with HASL called automaton-ABC. We apply this high-level method on several tasks of statistical inference and verification for biological CTMCs, including oscillatory models and time-bounded reachability problems. The implementation of our algorithms has led to a package in the Julia Programming language.

2. 1

 1 CTMC of a sleep cycle . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 1

 1 An illustration of the SMC ABC algorithm . . . . . . . . . . . . . . . 3.2 An illustration of Gaussian kernel density estimation. . . . . . . . . . 4.1 SIR Stochastic Petri Net with initial marking m 0 = (95, 5, 0) . . . . . 4.2 An example of LHA A count for SIR model. It accepts a trajectory if the number of infected people reaches eight. . . . . . . . . . . . . . . . 5.1 Plot of two 1-dimensional path. The L 1 distance of these paths is the sum of the areas hatched in green. . . . . . . . . . . . . . . . . . . . . 5.2 A trajectory from the simulated dataset. Left picture: 10-points observation. Right picture: 75-points observation. . . . . . . . . . . . . . 5.3 ABC-SMC posteriors with the ER model. N end = 10 6 . Left: timediscrete observations with 10 points, center: time-discrete observations with 75 points, right: the 5 trajectories. . . . . . . . . . . . . . . . . . 5.4 A trajectory from the simulated dataset. Left picture: 10-points observation. Right picture: 150-points observation. . . . . . . . . . . . . 5.5 ABC-SMC posterios with the SIR model. First row: N end = 10 6 , second row: N end = 10 7 . Left: time-discrete observations with 10 points, center: time-discrete observations with 75 points, right: the 5 trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Period automaton Aper. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Example of a oscillatory trajectory simulation synchronised with Aper. 5.8 Automaton-ABC posterior of r A with Aper automaton for the 1D experiment of doping 3-way oscillator. . . . . . . . . . . . . . . . . . . . 5.9 Correlation plot of automaton-ABC posterior with Aper automaton for the 3D experiment of doping 3-way oscillator. . . . . . . . . . . . . 5.10 Repressilator topology. Each protein represses the transcription of the successor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 5 simulated trajectories of the species P 1 . α varies in {50, 200, 1000, 4000}.102

  Aper posterior for the 3D experiment of repressilator model. . . . . . . . . . . . . . . . . 5.16 Correlation plot of automaton-ABC posterior with Aper posterior for the 4D experiment of repressilator model. . . . . . . . . . . . . . . . . 5.17 Automaton A ABC,ϵ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 40-points observation of the species P 1 from repressilator model over [0.0, 400.0] with a time step of 10.0. . . . . . . . . . . . . . . . . . . . . xv 5.19 Contour plot of ABC posteriors with marginals. Left: classical ABC posterior. Right: automaton-ABC posterior with A ABC,ϵ . . . . . . . . 6.1 Example of a simulation. In purple: the trajectory. In blue and green: the distances d 1 and d 2 . The rectangle is the F region. . . . . . . . . . 6.2 Examples of trajectories with zero-distance (top) and positive distance (bottom) from an F, a G and a U region (positive distances are depicted in red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Automaton for eventual property A F . . . . . . . . . . . . . . . . . . . 6.4 Automaton for global property A G . . . . . . . . . . . . . . . . . . . . . 6.5 Automaton for global and eventual property A G∧F when t 2 ≤ t 3 . . . . 6.6 Histogram of the automaton-ABC posterior with N = 500 particles. In red: estimated satisfaction probability function. In blue: the analytical satisfaction probability function. . . . . . . . . . . . . . . . . . . . . . 6.7 Trajectories of the ER system with θ left = (1, 1, 1) and time scale [0, 6], θ right = (0.1, 1, 0.1) and time scale [0, 30]. . . . . . . . . . . . . . . . . . 6.8 Trajectories of ER system for product P , with TR1, TR2 and TR3 regions and k 3 ∈ {10, 20, 50}. . . . . . . . . . . . . . . . . . . . . . . . 6.9 Average distances of the automata for the six formulae φ i , i ∈ {1, . . . , 6}, computed by the Statistical Model Checker Cosmos with approximation of 0.1 and 99% level of confidence. First row: A φ 1 , A φ 2 and A φ 3 . Second row: A φ 4 , A φ 5 and A φ 6 . . . . . . . . . . . . . . . . . . . . . . 6.10 Weighted histogram of automaton-ABC posteriors with 1000 particles. In each experiment, k 1 = k 2 = 1 and π k 3 (.) ∼ U (0, 100). In blue: the satisfaction probability function estimated on a selection of points with Prism model checker by numerical method. In red: the satisfaction function estimated through kernel density estimation method based on automaton-ABC samples. . . . . . . . . . . . . . . . . . . . . . . . xvi 6.11 Results of 2D experiments of ER system with 1000 particles (π k 1 (.), π k 2 (.) ∼ U (0, 100), k 3 = 1). Top: the 2D weighted histograms of the automaton-ABC posteriors. Middle: estimation of the satisfaction probability function with Prism by Statistical model checking (99% confidence interval with approximation 0.01). Bottom: kernel density estimation of the satisfaction probability function. . . . . . . . . . . . . . . . . . . . 6.12 Statistical Model Checking over [0.0, 0.0005] × [0.0, 20.0] with 10 points for the first axis and 20 points for the second. . . . . . . . . . . . . . .

N

  On the top left figure: weighted histogram of automaton-ABC posterior with 1000 particles for the 1D experiment; in blue: the true satisfaction probability function computed with the Prism model checker using the numerical engine of Prism over 40 points; in red: the estimated satisfaction function with the kernel density estimation method. The three other figures correspond to the 2D experiment. On the top right figure: the 2D histogram of the automaton-ABC posterior. On the bottom left figure: the kernel density estimation of the satisfaction probability function. On the bottom right figure: the estimation of the satisfaction probability function by the Prism model checker (numerical method). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 Results for the intracellular viral infection model with φ = G [0,50] G ≤ 10 ∧ F [50,200] G > 100. Left figure: the 2D histogram of the automaton-ABC posterior. Middle figure: kernel density estimation of the satisfaction probability function. Right figure: estimation of the satisfaction probability function by Monte-Carlo simulations (based on a 99% confidence level and approximation 0.01). . . . . . . . . . . . . . . . . 6.15 Estimation of the satisfaction probability functions in experiments TR1, TR2, TR3 for the ER system reported in Figure 6.10 with Smoothed MC algorithm. In blue: the estimated function; in green: the lower bound; in orange: the upper bound. . . . . . . . . . . . . . . A.1 Simulation of the SIR model with observation function g = [: I]. . . . A.2 Simulation of the ER model synchronised with a A F automaton. . . . A.3 Simulation of the ER model synchronised with a A F automaton. . . . A.4 UML graph of the different types defined in the package. . . . . . . . . xvii C.1 Weighted histograms of ABC. In green the true posterior distribution, in blue the true ABC posterior and in red the estimated ABC posterior with gaussian kernel. On the left: 1000 particles. On the right: 10000 particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Weighted histograms of ABC run. In green the true posterior distribution, in blue the true ABC posterior and in red the estimated ABC posterior with gaussian kernel. On the left: 1000 particles. On the right: 10000 particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 Histogram of ABC run with 1000 particles. In green the true posterior distribution, in blue the true ABC posterior and in red the estimated ABC posterior with gaussian kernel. On the left: multinomial resampling with the weights. On the right: weighted estimator. . . . . . . . of the ABC posterior for the three observations settings for the ER model with N = 10 6 particles. True value is k 3 = 5.0. . . . . . 5.2 Statistics of the ABC posterior for the three observations configurations for the SIR model with N = 10 6 particles. The true parameter is (k i , k r ) = (1.2E-3, 5E-2). . . . . . . . . . . . . . . . . . . . . . . . 5.3 Statistics of the ABC posterior for the three observations configurations for the SIR model with N = 10 7 particles. The true parameter is (k i , k r ) = (1.2E-3, 5E-2). . . . . . . . . . . . . . . . . . . . . . . . 5.4 Execution times of statistical inference with classical ABC compared to automaton-ABC with A ABC,ϵ . . . . . . . . . . . . . . . . . . . . . 6.1 Performance results for the one-dimensional experiments of automaton-ABC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Performance results for the two-dimensional experiments of automaton-ABC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Number of simulations before termination for both algorithms. For automaton-ABC: number of N = 1000 particles. For smoothed MC: each point of the dataset is estimated with 600 trajectories (default value). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Cosmos: Benchmark 1; 100 runs. . . . . . . . . . . . . . . . . . . . . . A.2 Cosmos: Benchmark 2; 5 runs. . . . . . . . . . . . . . . . . . . . . . . A.3 Cosmos: Benchmark 2 with 8 jobs (parallel execution); 5 runs. . . . . A.4 Catalyst.jl: Benchmark 1 . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Catalyst.jl: Benchmark 2 . . . . . . . . . . . . . . . . . . . . . . . . . xx List of Scripts1 Script (Simulation of the SIR model with SSA, Tau-leaping and ODE.) 2 Script (Bounded Kernel Density Estimation in Julia language) . . . . 3 Script (An example of a Cosmos run) . . . . . . . . . . . . . . . . . . 4 Script (Examples of ABC inference on parametric CTMC) . . . . . . . 5 Script (Oscillatory trends of parametric CTMC) . . . . . . . . . . . . 6 Script (Example CTMC) . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Script (A F LHA in Cosmos) . . . . . . . . . . . . . . . . . . . . . . . . 8 Script (A G LHA in Cosmos) . . . . . . . . . . . . . . . . . . . . . . . . 9 Script (A G∧F LHA in Cosmos) . . . . . . . . . . . . . . . . . . . . . . 10 Script (Automaton-ABC experiments related scripts) . . . . . . . . . . 11 Script (MarkovProcesses.jl package) . . . . . . . . . . . . . . . . . . . 12 Script (Simulation of the SIR model with MarkovProcesses.jl) . . . . . 13 Script (Synchronised simulation of the ER model with MarkovProcesses.jl) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Script (Experiment R1 of automaton-ABC with MarkovProcesses.jl) . 15 Script (Benchmarks of MarkovProcesses.jl) . . . . . . . . . . . . . . . xxi List of Listings A.1 Simulation of the SIR model. . . . . . . . . . . . . . . . . . . . . . . . A.2 Output of the print of the SIR variable created in Code A.1. . . . . . A.3 Creation of the SIR model with @network_model. . . . . . . . . . . . A.4 Output of the model created by a user in Code A.3. . . . . . . . . . . A.5 Simulation of the ER model synchronised with a A F automaton. . . . A.6 Automaton-ABC: Experience R1 of ER model with A F automaton. . (µ, σ 2 ) Normal distribution with mean µ and variance σ 2 P(λ) Poisson distribution with mean λ Exp(λ) Exponential distribution with mean 1 λ

  At the time of writing, the world has been hit for about a year by a pandemic of a new coronavirus disease called Covid-19. The pandemic has impacted society a lot,

  if a set of logical properties are satisfied by the stochastic process. The common representation of CTMCs used in this field is based on the transition rate matrix (see Remark 2.1.2). They are represented by a tuple M = (S, α, Q) where: • S is the set of states • α : S → [0, 1] with s∈S α(s) = 1 is the initial distribution. It corresponds to the distribution of S 0 , P S 0 .

Figure 2 .

 2 Figure 2.1 depicts an example of a CTMC with two states (S = {Sleep, Awake}).

Figure 2 . 2 :

 22 Figure 2.2: Graph of the SIR CTMC with the initial state s 0 = (95, 5, 0).

Figure 2 . 3 :

 23 Figure2.3: Effect of the population size on the simulations of the SIR model (only the number of infected people is plotted) with θ = (0.12/N pop , 0.05) and τ = 5.0. On the left plot, s 0 = (95, 5, 0) (N pop = 100). On the right plot, s 0 = (9500, 500, 0) (N pop = 10000). Red paths are simulated from the SSA. Green paths are simulated from the Tau-leaping algorithm. 5 paths for each algorithm.

Figure 2 . 4 :

 24 Figure2.4: Effect of the tau-leap time step on the simulations of the SIR model (only the number of infected people is plotted) with θ = (0.12/N pop , 0.05) and s 0 = (95, 5, 0) (N pop = 100). In the left plot, τ = 0.5. In the right plot, τ = 5.0. Red paths are simulated from the SSA. Green paths are simulated from the Tau-leaping algorithm.5 paths for each algorithm.

  ) i∈1:N are drawn from a multinomial distribution M with parameters N, (w m (x ) 1:N N samples with associated weights (w

Algorithm 12

 12 Simple ABC rejection Require: p(.) prior, ρ, η, ϵ Ensure: (θ (i) , y (i) ) 1≤i≤n drawn from the ABC posterior π ϵ ABC for i = 1 : n do repeat θ ∼ p(.) y ∼ p(.|θ) until ρ(η(y), η(y exp )) ≤ ϵ θ (i) , y (i) ← θ, y end for Remark 3.3.2 (ABC as a Bayesian inference with auxiliary likelihood) ABC can be seen as a standard Bayesian inference method with the likelihood p ABC (y exp |θ) = Y K ϵ (ρ(η(y), η(y exp )))p(y|θ)dy. An example based on Gaussian law is detailed in Appendix C. The approximate likelihood is analytically computed, which allows showing ABC in practice and test our algorithms.

  Figure 3.1: An illustration of the SMC ABC algorithm

  is moved locally through a proposal density qm (for example, a Gaussian kernel), as depicted in Figure 3.1. If the resulting simulation y ′ of θ (i) m verifies the acceptance condition ρ(η(y ′ ), η(y exp )), θ (i)

  Figure 3.2 illustrates this idea with a Gaussian kernel K.

  Figure 3.2: An illustration of Gaussian kernel density estimation.

  CTMC model checking suffer. These statistical methods are based on model simulations. The most classical one, called statistical model checking, is based on Monte Carlo simulations. Other statistical methods exist, and we propose a new one for parametric CTMCs in Chapter 6. More precisely, the goals of the Chapter are: • To describe the model checking formalism of a CTMC given a specification φ. • To define model checking of parametric CTMC, which is model checking on a collection of parametrised CTMCs.

•

  To give an overview of statistical methods related to model checking.• To detail the HASL formalism, which is well-suited for CTMC verification.The Chapter is organised as follows. Section 4.1 describes temporal logics to formally define model specifications. Section 4.2 describes the verification tasks for a CTMC and the related statistical methods. Section 4.3 describes the verification tasks for a collection of CTMCs and the related statistical methods. Section 4.4 describes the HASL formalism based on Linear Hybrid Automata.

  MITL formulae do not involve any probabilistic operator. The truth of an MITL formula for a path σ ∈ P ath(M) is expressed through a satisfaction relationship denoted |= and defined below (see Section 2.2.1 of Chapter 2 for the notations of a CTMC path).

  Definition 4.1.2 (Satisfaction relation |= of MITL)

  Definition 4.1.4 (Satisfaction relation |= for a CSL state formula ϕ)

  Satisfaction probability of a path formula) Let M be a CTMC and φ be an MITL formula. The probability that M satisfies φ defined by P r(φ; M) = P r M ({σ ∈ P ath(M), σ |= φ}) is called the satisfaction probability of φ. Two kinds of problems related to model checking are distinguishable. Definition 4.2.2 (Estimation problem in model checking for CTMC) Given an MITL formula φ and a CTMC M, the estimation problem in model checking for CTMC is defined as: Compute P r(φ; M) Definition 4.2.3 (Threshold problem in model checking for CTMC) Given a CSL formula ϕ = P ▷◁p φ with φ a path formula, and a CTMC M, the threshold problem in model checking for CTMC is defined as: Verify M |= ϕ, i.e. P r(φ; M) ▷◁ p. Remark 4.2.1 (Qualitative/quantitative approach of Probabilistic Model Checking)

Chapter 4 .

 4 Verification of Continuous-Time Markov Chains 67 With this equality, one can set the number N of simulations according to the radius r = -log(α/2) 2N

  a vector θ naturally parametrises CTMCs. CTMCs induced by Chemical Reaction Networks (Section 2.3) are an illustrative example: they are parametrised by reaction rate constants. In this section, we present verification methods for a collection of CTMCs indexed by a parameter. These collections are called parametric Continuous-Time Markov Chains, for which we give a definition below. Definition 4.3.1 (Parametric CTMC)

Figure 4 .Figure 4 . 1 :

 441 Figure 4.1 illustrates the SIR CRN with the SPN formalism.

Figure 4 . 2 :

 42 Figure 4.2: An example of LHA A count for SIR model. It accepts a trajectory if the number of infected people reaches eight.

  y) | M ax(y) | Int(y) | M ean(y) y ::= c | x | y + y | y • y | y/y In this definition, y is an arithmetic expression built on top of the automaton variables, computed at each step. Y is called an HASL trajectory expression. It represents a value computed at the end of an accepted trajectory. Z is an HASL expression, evaluated over a set of simulated trajectories. For Z expressions, P ROB is the probability that a path is accepted by the LHA, AV G(Y ) is the first moment (expectation) of the HASL trajectory expression Y , which allows the computation of any k-th order moments. For an expression Y , several quantities are available: • Last(y) is the value of y in the last state • M in(y), M ax(y), M ean(y) are the minimum, maximum and average value of y along the trajectory • Int(y) is the integral over time of y along the trajectory Definition 4.4.4 (Notation for synchronised simulations) Let Y an HASL trajectory expression, A an LHA and M a CTMC. The notation x ∼ (Y, A) × M describes the evaluation of the trajectory expression Y after a simulation of M synchronised with A. Considering the synchronised simulation example of Section 4.4.2, x ∼ (Last(n), A count )× M means that x stores the last value of n, which is 8. The notation of Definition 4.4.4 is mainly used in the automaton-ABC procedure detailed in the next chapter.

  Statistical model of a parametric CTMC) Let (S θ ) θ∈Θ a parametric CTMC. Then (P S θ ) θ∈Θ is a statistical model dominated by a measure µ. Corollary 5.1.1 (Density of a CTMC)

  dz 0:K Chapter 5. Automaton-ABC for the statistical inference of CTMCs 85 When the true state is observed at each step, Y k = S t k and the likelihood becomes l

  , they performed exact posterior distribution inference for the repressilator model (a model considered in the last two Sections of the Chapter) using the Sequential Importance Sampling algorithm. It lasted 7137 hours after a state-space reduction. The author had difficulties managing the run of the Particle Marginal 5.1. Observation model and likelihood Metropolis-Hasting algorithm, because assessing a trustful convergence of a Markov Chain is a challenging task. Driving the convergence of ABC-SMC (Algorithm 14) was easier because the Sequential Monte Carlo scheme adapts the tolerance ϵ during the run, and the α adaptive tolerance parameter (Section 3.3.2) easily adapts the convergence speed.

Figure 5 . 1 :

 51 Figure 5.1: Plot of two 1-dimensional path. The L 1 distance of these paths is the sum of the areas hatched in green.

Figure 5

 5 Figure 5.2: A trajectory from the simulated dataset. Left picture: 10-points observation. Right picture: 75-points observation.

Figure 5 .

 5 Figure 5.3 shows the resulting ABC-SMC posteriors.Table 5.1 reports statistics

Figure 5

 5 Figure 5.4: A trajectory from the simulated dataset. Left picture: 10-points observation. Right picture: 150-points observation.

Figure 5 .

 5 Figure5.5 shows the resulting ABC-SMC posteriors. Tables 5.2 and 5.3 report statistics about the posteriors. For the two simulation budgets N end = 10 6 and N end = 10 7 , no significant improvement in the parameter inference has been noticed for the 10-points observations than the two more informative configurations.

  3) with HASL formalism called automaton-ABC. As classical ABC for a parametric CTMC (M θ ) θ∈Θ , automaton-ABC requires a prior π and a tolerance ϵ. In this context, we do not consider observations y exp . Instead, we consider an HASL trajectory expression Y and an LHA A. The main difference with classical ABC is that the simulation of the CTMC is replaced by the synchronised simulation (Y, A) × M (Definition 4.4.4). Algorithm 15 details the automaton-ABC analogous to the simple ABC Algorithm 12. Chapter 5. Automaton-ABC for the statistical inference of CTMCs 93 Algorithm 15 General automaton-ABC Algorithm Require: (M θ ) θ∈Θ parametric CTMC, π prior, N : number of particles, ϵ: tolerance level, Y : HASL trajectory expression, A : LHA Ensure: (θ

Figure 5 . 7 :

 57 Figure 5.7: Example of a oscillatory trajectory simulation synchronised with Aper.

Figure 5 .

 5 Figure 5.7 illustrates the behaviour of Aper for a trajectory simulated by the doping 3-way oscillator (see section below) with L = 300, H = 360, N per = 5. Let us describe this trajectory step by step:

Figure 5 .

 5 Figure 5.8 shows the resulting automaton-ABC posterior. The sequential run lasted 140 seconds and performed 19846 simulations. First, we can observe that the support is included in [0.0, 4.0], which is shrinker than the support of the prior:

Figure 5 . 8 :

 58 Figure 5.8: Automaton-ABC posterior of r A with Aper automaton for the 1D experiment of doping 3-way oscillator.

Figure 5 . 9 :

 59 Figure 5.9: Correlation plot of automaton-ABC posterior with Aper automaton for the 3D experiment of doping 3-way oscillator.

Figure 5

 5 Figure 5.10: Repressilator topology. Each protein represses the transcription of the successor.

Figure 5 .

 5 Figure 5.11 and Figure 5.12 show simulations with different parameter values of α and β. It illustrates trajectories with different amplitudes and period duration, but oscillations are pretty stable.

Figure 5 .

 5 Figure 5.11: 5 simulated trajectories of the species P 1 . α varies in {50, 200, 1000, 4000}.

Figure 5 .

 5 Figure 5.12: 5 simulated trajectories of the species P 1 . β varies in {0.5, 1.0, 2.0, 4.0}.

Figure 5 .

 5 Figure 5.13 shows simulations with different parameter values of n. This parameter has more influence on the stability of oscillations since simulations with n = 0.5 and n = 1.0 are very noisy.

Figure 5 .

 5 Figure 5.13: 5 simulated trajectories of the species P 1 . n varies in {0.5, 1.0, 2.0, 5.0}.

Figure 5 .

 5 Figure 5.14: 5 simulated trajectories of the species P 1 . α 0 varies in {0.0, 0.01, 0.1, 1.0}.

Figure 5 .

 5 Figure 5.15 shows the correlation plot of the resulting automaton-ABC posterior. The parallel run (250 jobs) performed 15333 simulations and lasted 2040 seconds. The parameter n has the shrinkest marginal posterior, consistent with the simulation study detailed above: varying n induces more instability than α and β.

Figure 5 . 15 :

 515 Figure 5.15: Correlation plot of automaton-ABC posterior with Aper posterior for the 3D experiment of repressilator model.

Figure 5 . 16 :

 516 Figure 5.16: Correlation plot of automaton-ABC posterior with Aper posterior for the 4D experiment of repressilator model.

Figure 5 .

 5 Figure 5.18: 40-points observation of the species P 1 from repressilator model over [0.0, 400.0] with a time step of 10.0.

Figure 5 .

 5 Figure 5.19 shows the ABC posteriors for both ABC procedures. These two posteriors are very similar, which means that the two procedures behaved in the same way, as expected.

Figure 5 . 19 :

 519 Figure 5.19: Contour plot of ABC posteriors with marginals. Left: classical ABC posterior. Right: automaton-ABC posterior with A ABC,ϵ .

  (ABC) methods (Section 3.3) allow approximating the posterior distribution of a model without evaluating the likelihood function when the computation cost is too high or even impossible. These methods rely on model simulations and a discrepancy between simulations and observations. Simply speaking, only parameters for which simulated summary statistics are close to observed ones, relative to a distance, are preserved while the others are dismissed.These parameters are sampled from the ABC posterior, which approximates the true posterior distribution. This chapter proposes a new method based on the ABC algorithm for the estimation problem of parametric CTMC in verification, described in Section 4.3.1. We want to estimate the satisfaction probability function given a time-bounded reachability specification φ related to a parameter θ. The initial idea of this work relies on transposing this concept originally developed with observations data to temporal logic: only parameters for which simulations are close to fulfilling a given formula are kept while the others are discarded.

  t≥t 2 ∧γ,U ------→ defines two edges:♯,t≥t 2 ∧γ,U -------→ and ♯,E(x)=0∧γ,U

Figure 6 . 3 :

 63 Figure 6.3: Automaton for eventual property A F .

Figure 6 . 4 :

 64 Figure 6.4: Automaton for global property A G .

Figure 6 . 5 :

 65 Figure 6.5: Automaton for global and eventual property A G∧F when t 2 ≤ t 3 .

  are drawn from the prior π, and the first tolerance level ϵ to reach equals the α-quantile of the distances d i , resulting from the synchronised simulations M θ (i) 0 × A φ . Then, at each iteration m, a parameter is drawn among the parameters θ (i) m (i ∈ {1, . . . , N }), and is moved by a kernel distribution K, and the synchronised simulation M θ (i) m × A φ is performed. This procedure is done until the resulting distance d i is below the current tolerance level ϵ, inducing that the parameter θ (i) m is kept. After doing so for the N parameters, we compute a new tolerance level ϵ that equals the α-quantile of the distances d i . These iterations are repeated until the last tolerance level ϵ = 0 is reached. The introduction of several steps with positive decreasing tolerances leads to an efficient exploration of the parameter space driven by A φ . Another interesting point is that any algorithm inherited from the family of ABC methods keeps the guarantee of Theorem 6.4.1. Thus depending on the problem, one can use any Monte-Carlo based sampler (ABC-MCMC, ABC-SMC and so on) and still estimate the satisfaction probability function.

Figure 6 . 6 :

 66 Figure 6.6: Histogram of the automaton-ABC posterior with N = 500 particles. In red: estimated satisfaction probability function. In blue: the analytical satisfaction probability function.

Figure 6 .

 6 Figure 6.7 shows two (4-dimensional) trajectories sampled from the CTMC model M θ of the enzymatic reaction system with parameters θ = (1, 1, 1) (left) and θ = (0.1, 1, 0.1) (right).

Figure 6 . 7 :

 67 Figure 6.7: Trajectories of the ER system with θ left = (1, 1, 1) and time scale [0, 6], θ right = (0.1, 1, 0.1) and time scale [0, 30].

Figure 6 . 8 :

 68 Figure 6.8: Trajectories of ER system for product P , with TR1, TR2 and TR3 regions and k 3 ∈ {10, 20, 50}.

Figure 6 .

 6 Figure 6.8 shows batches of trajectories of the ER model for product P , corresponding to parameters θ 1 : (1, 1, 50), θ 2 = (1, 1, 20) and θ 3 = (1, 1, 10). It appears that trajectories for M θ 1 (red) are more likely to traverse T R 1 , those for M θ 2 (blue) to traverse T R2 and those for M θ 3 (green) to traverse T R3.

Figure 6 . 9 :

 69 Figure 6.9: Average distances of the automata for the six formulae φ i , i ∈ {1, . . . , 6}, computed by the Statistical Model Checker Cosmos with approximation of 0.1 and 99% level of confidence. First row: A φ1 , A φ2 and A φ3 . Second row: A φ4 , A φ5 and A φ6 .

Figure 6 . 10 :

 610 Figure 6.10: Weighted histogram of automaton-ABC posteriors with 1000 particles. In each experiment, k 1 = k 2 = 1 and π k3 (.) ∼ U (0, 100). In blue: the satisfaction probability function estimated on a selection of points with Prism model checker by numerical method. In red: the satisfaction function estimated through kernel density estimation method based on automaton-ABC samples.

Figure 6 .

 6 Figure 6.10 illustrates the evaluation of the posterior distribution π φ-ABC and the estimation of the satisfaction probability as a function of parameter k 3 (1D experiments) obtained with the automaton-ABC method (Algorithm 18) over the formulae φ 1 , φ 2 , φ 3 .

Figure 6 .F

 6 Figure 6.11 depicts the results of the 2D experiments on the ER system with φ 4 , φ 5 , φ 6 formulae. The triangular profile of the joint posterior in experiments T R4 and T R5 (computed with k 3 = 1 and π k 1 (.), π k 2 (.) ∼ U (0, 100)) indicates that only very low values of k 1 (k 1 ≤ 0.015 for T R4, k 1 ≤ 1 for T R5) combined with rather high

Figure 6 .

 6 Figure 6.12: Statistical Model Checking over [0.0, 0.0005]×[0.0, 20.0]with 10 points for the first axis and 20 points for the second.

Figure 6 . 13 :

 613 Figure 6.13: Results for the SIR model with φ = G [0,100] (I > 0) ∧ F [100,120] (I = 0). On the top left figure: weighted histogram of automaton-ABC posterior with 1000 particles for the 1D experiment; in blue: the true satisfaction probability function computed with the Prism model checker using the numerical engine of Prism over 40 points; in red: the estimated satisfaction function with the kernel density estimation method. The three other figures correspond to the 2D experiment. On the top right figure: the 2D histogram of the automaton-ABC posterior. On the bottom left figure: the kernel density estimation of the satisfaction probability function. On the bottom right figure: the estimation of the satisfaction probability function by the Prism model checker (numerical method).

Figure 6 . 14 :

 614 Figure 6.14: Results for the intracellular viral infection model with φ = G [0,50] G ≤ 10 ∧ F [50,200] G > 100. Left figure: the 2D histogram of the automaton-ABC posterior. Middle figure: kernel density estimation of the satisfaction probability function. Right figure: estimation of the satisfaction probability function by Monte-Carlo simulations (based on a 99% confidence level and approximation 0.01).

  This computational time has to be put in perspective to classical Statistical Model Checking methods. For example, in experiment T R4, Statistical Model Checking has been performed over a grid of 200 points of a much smaller set than [0, 100] × 6.6. A comparison with Smoothed Model Checking reproduced the 1D experiments for the ER model (Section 6.5.2) using the Smoothed MC tool.

Figure 6 .

 6 Figure6.15 depicts plots of the satisfaction probabilities obtained with Smoothed MC, which show a good agreement with those obtained with the automaton-ABC approach and Prism model checker (Figure6.10). Table6.3 reports the number of trajectories generated by both algorithms, showing a clear advantage for the automaton-ABC approach. In this table, we omit simulations involved in estimating P r(φ; M θ * ) for the automaton-ABC. On the contrary, this table does not show the cost of normalisation constant estimation in the Smoothed MC posterior (Equation4.1) by Expectation-Propagation. Also, a default value of 600 trajectories is set for each satisfiability probability estimation by Statistical Model Checking, which should be higher if one wants a high confidence level and a precise approximation.

Figure 6 . 15 :

 615 Figure 6.15: Estimation of the satisfaction probability functions in experiments TR1, TR2, TR3 for the ER system reported in Figure 6.10 with Smoothed MC algorithm. In blue: the estimated function; in green: the lower bound; in orange: the upper bound.

  Another solution is to defined the set R O ⊆ ALL of reactions that change the number of the studied species O. If R ∈ R O , the species O is a reactant or a product of the reaction R (its population number changes). If R ∈ ALL \ R O occurs, the number of species O remains equal. Then, we only have to change ALL by R O in the automata edges. This also improves the computational cost, because n and the distance d are updated only when a reaction in R O occurs. model a wide range of biological phenomena. If traditionally their dynamics are described by Ordinary Differential Equations, this approach is not relevant when the population is low. It is especially the case when one studies interactions between molecular species at the scale of the living cell (genetic regulatory networks, signalling pathways). By the reduction of such assumption, stochastic dynamics of CRN are described by CTMCs.

7. 1 .

 1 Limits and Perspectives specification. In this context, it was used to stop the simulation each time the Euclidean distance to observations exceeds the tolerance ϵ in the ABC procedure. In an experiment with the repressilator model, we showed that ABC computational time was divided by more than three compared to classical ABC. Fourth, we used automaton-ABC for time-bounded reachability problems. The satisfiability distances represent the distance of a trajectory σ from specific forms of MITL formula φ. We designed LHA to compute these distances and used automaton-ABC to address quantitative model checking for parametric CTMCs. More precisely, we highlighted the link between the satisfaction probability function P r(φ; M) and the automaton-ABC posterior with A φ named π φ-ABC . Automaton-ABC-SMC allows exploring the parameter space concerning the formula φ efficiently. The results were obtained over several models of CRN and constantly compared to either an analytical solution, model checking with Prism, or statistical model checking with Cosmos or Prism. Due to optimisation and parallelisation of our implementation, the runs of automaton-ABC did not exceed two minutes, and outperformed state-of-the-art statistical model checking methods.

Figure A. 2 :

 2 Figure A.2: Simulation of the ER model synchronised with a A F automaton.

Figure A. 2

 2 Figure A.2 can be generated by running the Julia script located in code/appendix_ pkg/sim_sync_er.jl.

Figure C. 2 :

 2 Figure C.2: Weighted histograms of ABC run. In green the true posterior distribution, in blue the true ABC posterior and in red the estimated ABC posterior with gaussian kernel. On the left: 1000 particles. On the right: 10000 particles.

Titre:

  Inférence et vérification statistiques de réseaux de réactions chimiques. Mots clés: Réseaux de réactions chimiques, inférence, vérification, ABC, HASL, Julia Résumé: Les réseaux de réactions chimiques (CRN) constituent un formalisme utilisé pour modéliser des processus biologiques. Sous certaines hypothèses, le processus stochastique sous-jacent pour décrire ses dynamiques est une chaîne de Markov en temps continu (CTMC): l'état futur du système ne dépend que de l'état courant. L'inférence statistique de ce type de CTMC est complexe: le calcul de la vraisemblance est en général difficile à résoudre. Les méthodes ABC (Approximate Bayesian Computation) forment une classe de méthodes bayésiennes sans calcul de vraisemblance qui permettent d'approcher la distribution postérieure avec des simulations de Monte-Carlo. La vérification de modèles se penche de plus en plus sur la biologie des systèmes. En effet, il y a un réel besoin de comprendre les interactions complexes entre molécules dans les systèmes biologiques. Des méthodes de vérification statistiques ont été développées pour appliquer ces méthodes dans des modèles plus complèxes. Le principe est de simuler un certain nombre de fois le modèle et de calculer le ratio des simulations qui ont vérifié une propriété. Récemment, une logique temporelle appelée HASL a été introduite pour la vérification statistique de modèles: elle adopte intrinsèquement le point de vue statistique de la vérification. Dans cette thèse, nous nous intéressons à l'inférence statistique et la vérification statistique de chaines de Markov en temps continu définies par un modèle de réseaux de réactions chimiques. Notre contribution tient principalement dans la formulation d'un algorithme ABC combiné avec le formalisme HASL appelé automaton-ABC. Nous appliquons cette méthode haut niveau sur plusieurs tâches d'inférence statistique et de vérification pour des CTMCs issus de systèmes biologiques, impliquant notamment des modèles oscillatoires et des problèmes d'atteignabilité bornés en temps. L'implémentation des méthodes présentées est rendue disponible sous la forme d'une bibliothèque dans le langage de programmation Julia. Title: Statistical inference and verification of Chemical Reaction Networks Keywords: Chemical Reaction Networks, inference, model checking, ABC, HASL, Julia

Chapter 2 Markov Chains and Chemical Reaction Networks

This chapter is an introduction to Markov Chains theory. The idea is to present Markov Chains with stochastic process and model checking point of views and see how they describe the dynamics of Chemical Reaction Networks, a formalism that defines the interactions between chemical species. Markov Chains are a specific class of stochastic models. They are characterised by the memoryless property: if a value is known at a specific time, the future is independent of the older values. In our work, this class of models is involved in several tasks: • Continuous-Time Markov Chains, also called Markov Jump Processes, is the underlying probabilistic model involved in Chemical Reaction Networks, our main study subject. • Discrete-Time Markov Chains, sometimes simply called Markov Chains, are involved in advanced statistical methods (Markov Chain Monte Carlo, Hidden Markov Models).

  Let (S n ) n∈N be an ergodic Markov Chain with n-step transition probabilities p

	Theorem 2.1.2 (Ergodic theorem)
	(n) ij and
	its stationary distribution π. Then the n-step transition probabilities converge to the
	stationary distribution, i.e
	∀i, sup

j

  1 16 . In this CTMC, only one edge arises from the Sleep state. So the sojourn time of the Sleep state is T Sleep ∼ Exp(Q(Sleep, Awake)). With a time scale of hours, this modelisation considers that the average time of sleep is 8 hours.

	1
	16
	Sleep	Awake
	1
	8
	Figure 2.1: CTMC of a sleep cycle
	Depending on the literature, the property s ′ ∈S Q(s, s ′ ) = 0 can be either in-
	cluded in the definition or not. We assume this property to hold. Simply speaking,
	we do not allow self transitions (s

t

-→ s).

  y|y exp )dy Chapter 3. Statistical methods 47 It shows that the smaller the tolerance level ϵ is, the closer we are from the posterior. It is then essential to understand the different levels of approximations of

	ABC:
	π ϵ ABC (.|η(y exp )) → (1)	π ϵ ABC (.|η(y exp )) → (2)	p(.|η(y exp )) ≈ (3)	p(.|y exp )
	1. How the empirical posterior π ϵ ABC (.|η(y exp )) based on the samples (θ i , y i ) i is
	close to π ϵ ABC ?
	2. How to sample efficiently with small tolerance ϵ to be close to p(.|η(y exp ))?
	3. How informative should be η for p(.|η(y exp )) to be a good approximation of the
	true posterior p(.|y exp )?
	Remark 3.3.1 (Proposal density and kernel function)
	lim ϵ→0	π ϵ ABC (θ|y exp ) ∝ lim ϵ→0 Y	π ϵ ABC (θ, y|y exp )dy
			∝ lim ϵ→0 Y	K ϵ (ρ(y, y exp ))p(y|θ)p(θ)dy
			∝

Y δ yexp (y)p(y|θ)p(θ)dy ∝ p(y exp |θ)p(θ)

  δ s , Q, AP, L)

	s |= S ▷◁p ϕ	⇔	lim

t→+∞ P r M ({σ ∈ P ath(M), σ@t |= ϕ}) ▷◁ p P r M 0 ({σ ∈ P ath(M 0 ), σ |= ϕ}) is the probability that σ |= ϕ when the CTMC M begins at state s. Remark 4.1.1 As detailed in the Remark 2.2.2 of Chapter 2, P r M is often denoted P r α in the literature. With this notation, P r M 0 ({σ ∈ P ath(M 0 ), σ |= ϕ}) = P r s ({σ ∈ P ath(s), σ |= ϕ}). Definition 4.1.5 (Satisfaction relation |= for a CSL path formula φ)

  have a good estimation over Θ, observations on the dataset are added iteratively (i.e. drawing a new θ * and perform statistical model checking over M θ * ) until a convergence criterion is met (about the estimated error).As the matter of fact, Smoothed MC is not Gaussian Process regression. In classical GP regression, we consider observations (x i , y i ) 1:N where y i ∼ Y x i and Y x 1:N = (Y x 1 , . . . , Y x N ) t as the random sample. Let x * ∈ X . We want to address the prediction task, i.e. estimate f (x

	Remark 4.3.2 (Smoothed MC and Gaussian Process regression)

* ). This task is done by identifying the conditional law Y x * |Y x 1:N . As Y is a Gaussian process, we know Y x * , Y x 1 , . . . , Y x N is a Gaussian vector:

Table 5

 5 ||.|| 2 , 10 points d ||.|| 2 , 75 points d L 2

	d ||.|| 2 , 10 points	d ||.|| 2 , 75 points	d L 2
	Figure 5.3: ABC-SMC posteriors with the ER model. N end = 10 6 .
	Left: time-discrete observations with 10 points, center: time-discrete
	observations with 75 points, right: the 5 trajectories.

.1 reports statistics about the posteriors. These results indicate that with the same simulation budget, the resulting ABC posteriors are pretty similar in estimation accuracy. The estimated mean (Table

5

.1, first row) with the 10-points observations is slightly better, even if it contains less information than the two other settings.

Table 5 . 1 :

 51 Statistics of the ABC posterior for the three observations settings for the ER model with N = 10 6 particles.True value is k 3 = 5.0.

Table 5 .

 5 ||.|| 2 , 150 points d L 2 Table 5.3: Statistics of the ABC posterior for the three observations configurations for the SIR model with N = 10 7 particles. The true parameter is

	Mean k i	1.25E-3	1.27E-3	1.27E-3
	Std k i	3.4E-4	2.90E-4	3.25E-4
	Mean k r	4.57E-2	4.75E-2	4.73E-2
	Std k r	6.96E-3	6.54E-3	6.6E-3

2: Statistics of the ABC posterior for the three observations configurations for the SIR model with N = 10 6 particles. The true parameter is (k i , k r ) = (1.2E-3, 5E-2). d ||.|| 2 , 10 points d

Table 5 .

 5 4 shows the execution time of both runs. The automaton-ABC run was at least three times faster than the classical ABC procedure.

		Number of jobs Number of simulations Execution time (sec)
	Automaton-ABC	80	2300706	16294
	Classical ABC	80	2554773	54049

Table 5 .

 5 4: Execution times of statistical inference with classical ABC compared to automaton-ABC with A ABC,ϵ

  1 , t 2 ]. These distances are related to the time-bounded reachability problems described in Section 6.1. The following definitions hold for any trajectory σ, whose last jump occurred at time t last ≤ t 2 , wheret last = min{t ∈ R ≥0 , ∀t ′ ∈ [t, t 2 ], σ@t ′ = σ@t}. -s ′ [i]) 2 isthe least Euclidean distance between s and S µ , whereas d e ((t, s), S µ [t 1 ,t 2 ] ) is the least Euclidean distance between (t, s) and Let M a CTMC, [t 1 , t 2 ] ⊆ R ≥0 . Let σ ∈ P ath(M) a trajectory observed until t 2 , with

	Indeed, jumps after t 2 do not affect a reachability problem over [t 1 , t 2 ]. In the fol-
	lowing, d e (s, S µ ) = min s ′ ∈Sµ i=1 (s[i] S µ d [t 1 ,t 2 ] .
	Definition 6.2.4 (Distances from several MITL formulas)

  the time-bounded satisfiability regionS [t 1 ,t 2 ] (x 1 ≤x O ≤x 2 ) .It uses the same variables as A F plus an extra timer t ′ , to measure the duration of a segment falling outside the satisfiability region within [t 1 , t 2 ], and a boolean flag in, which is set to true if the last state of the path is originating in [t 1 , t 2 ] and is outside of the region [x 1 , x 2 ]. in is used to distinguish cases where the path is out of the region [x 1 , x 2 ] with t ′ < t 1 , and a new event occurs at time t ′′ > t 1 , to add

Remerciements

Code A.1 simulates the SIR model and plots Figure A.1. Let us describe this code line by line. i. using MarkovProcesses loads the package.

Code A.3 creates a SIR model and stores it in my_model . Each reaction is described by a couple of expressions. The first expression describes the interactions between species (reactants => products). The second expression describes the kinetic rate (Equation 2.4). Each symbol that cannot be found in any reaction is considered as a parameter. The output of println(my_model) is given in Listing A.4.

Summary

• We have applied the automaton-ABC procedure on two main tasks. First, we describe oscillatory trends of models such as repressilator with the help of A per .

Second, we have speeded up classical ABC inference with the help of A ABC,ϵ .

The next chapter tackles the estimation problem in model checking parametric CTMC (Section 4.3.1) with the help of automaton-ABC.

6.1. Problem setting: time-bounded reachability leading to an efficient exploration of the parameter space. Our method also allows the estimation of the satisfaction function of a formula, with a remarkable result linking this satisfaction function to the ABC posterior (Theorem 6.4.1). We demonstrate the effectiveness of our method by its application to several models of Chemical Reaction Networks (CRN). This Chapter is organised as follows. In Section 6.2, the notion of satisfiability distance for time-bounded reachability problems is introduced. We see in Section 6.3 how simulations synchronised with Linear Hybrid Automata compute the satisfiability distances. The automaton-ABC procedure with LHA related to satisfiability distances is detailed in Section 6.4: the objective is to find the parameter subspace for which the probability of reaching the target region is positive. Its effectiveness is demonstrated through several experiments in Section 6.5 and Section 6.6, while some conclusive remarks and perspectives are discussed in Section 6.7.

Problem setting: time-bounded reachability

In this Chapter, we tackle the estimation problem detailed in Section 4.3.1. We consider a parametric CTMC (M θ ) θ∈Θ (Section 4.3.1) with a state space of dimension d. This work is focused on CTMCs whose states represent the number of individuals of several species, which is also called Markov Population Processes (see Remark 2.3.1). φ is a specification described by an MITL formula (Definition 4.1.1). Our goal is to estimate the satisfaction probability function f φ (Definition 4.3.2):

where P r(φ; M θ ) (Definition 4.2.1) is the probability that a trajectory of M θ verifies φ.

In particular, we consider specifications φ related to time-bounded reachability.

The term reachability problem identifies the class of problems that checks if a given model reaches (i.e. enters) at some point during its execution a specific region of its state space. For models that produce time series, time-bounded reachability aims to establish whether the state space's target region is entered within a time-interval

Based on MITL formulae, we distinguish three kinds of time-bounded reachability problems: Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs

113

• Eventual reachability problems (F I µ) are concerned with model trajectories that fulfil µ at least once within a given time interval I.

• Global reachability problems (G I µ) are concerned with trajectories that consistently fulfil µ within a time interval I.

• Conditional reachability problems (µ 1 U I µ 2 ) are concerned with trajectories that fulfil µ 2 at least once within a time interval I and fulfil µ 1 beforehand.

Each kind of problem induces a satisfiability region that characterises the trajectories verifying the corresponding formula. A region is elementary if it is characterised by a single subset of S. We denote [[p, q]] the set of integers {p, . . . , q -1, q}. For example, let a bi-dimensional CTMC S with a state space S = N 2 . R 1 = [[1, 2]]×N is an elementary region (S[1] in [[1, 2]] while S[2] is unbounded), R 2 = ([[0, 3]] ∪ [[5, 8]]) × [[5, ∞[[ is a non-elementary region (S[1] is either in [[0, 3]] or [[5, 8]], S [2] is larger than 5), whereas T R 1 = ( [[1, 2]]×N) ×[0.2, 1.41] is an elementary time-bounded region (similar to R 1 , but with the supplemental condition that the time is in [0.2, 1.41]).

Satisfiability distances

The adaptation of the parameter space exploration by ABC algorithms to reachability problems requires the definition of a satisfiability distance d (σ, φ). This distance represents how far a trajectory σ ∈ P ath(M) is from satisfying an MITL formula φ. Note that we call it a distance to keep up with the ABC vocabulary: it is not mathematically strictly a distance function, but rather a discrepancy. Indeed, (σ, φ) is not a pair of elements of the same set.

Definition 6.2.1 (Satisfiability distance)

Let φ an MITL formula. f is a satisfiability distance if ∀σ, f (σ, φ) ≥ 0 and f (σ, φ) = 0 ⇔ σ |= φ 6.2. Satisfiability distances . Otherwise, the distance d(σ, F [t 1 ,t 2 ] µ) is the minimal Euclidean distance of any point of σ within [t 1 , t 2 ] from S µ . d(σ, G [t 1 ,t 2 ] µ) is the integral of the Euclidean distance from S µ of any point of σ within [t 1 , t 2 ]. The bigger the distance, the more the trajectory is globally outside of

Concerning the Until formula, the distance (6.3) bears three components:

• d(σ, G [0,t 1 ] µ 1 ) accounts for the fact that a trajectory satisfying (µ 1 U [t 1 ,t 2 ] µ 2 ) must never leave region µ 1 before t 1 .

• d(σ, F [t 1 ,t 2 ] µ 2 ) accounts for the fact that σ must enter region µ 2 within [t 1 , t 2 ]

• d(σ, G [t 1 ,t min ] (µ 1 ∨ µ 2 )) accounts for the fact that there must be a time t min where σ switches from region µ 1 to region µ 2 directly, without spending time in an intermediate region. If that is not the case (i.e. if σ within [t 1 , t 2 ] has points in the complementary region ¬(µ 1 ∨ µ 2 )), then the distance 6.3 is incremented by this term that quantifies how far we are from one of the two regions µ 1 or µ 2 . Remark 6.2.1 (Until distance)

In the following, we only designed an LHA (A G∧F ) for computing the Until distance when µ 1 ∨µ 2 (the exclusive disjunction) is always true over any state of the CTMC, which implies the third term of the distance is zero. Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs 117 distance = 0 t 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " Y R 9 G K J d m / n A V x Y n m S I + 1 k t a J 3 f 0 = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l q Q Y 8 F L x 4 r 2 g 9 o Q 9 l s N + 3 S z S b s T o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I J H C o O t + O 4 W N z a 3 t n e J u a W / / 4 P C o f H z S N n G q G W + x W M a 6 G 1 D D p V C 8 h Q I l 7 y a a 0 y i Q v B N M b u d + 5 4 l r I 2 L 1 i N O E + x E d K R E K R t F K D z j w B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p x B U y S Y 3 p e W 6 C f k Y 1 C i b 5 r N R P D U 8 o m 9 A R 7 1 m q a M S N n y 1 O n Z E L q w x J G G t b C s l C / T 2 R 0 c i Y a R T Y z o j i 2 K x 6 c / E / r 5 d i e O N n Q i U p c s W W i 8 J U E o z J / G 8 y F J o z l F N L K N P C 3 k r Y m G r K 0 K Z T s i F 4 q y + v k 3 a t 6 l 1 V a / f 1 S q O e x 1 G E M z i H S / D g G h p w B 0 1 o A Y M R P M M r v D n S e X H e n Y 9 l a 8 H J Z 0 7 h D 5 z P H w K w j Z I = < / l a t e x i t > t 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " q 4 Z 1 m j f / m 0 K X 4 d q 2 z k H w f E B D B K I = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l q Q Y 8 F L x 4 r 2 g 9 o Q 9 l s N + 3 S z S b s T o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I J H C o O t + O 4 W N z a 3 t n e J u a W / / 4 P

< l a t e x i t s h a 1 _ b a s e 6 4 = " y F e 8 s e V d C S C u p 3 8 I r / J y h 1 j A

< l a t e x i t s h a 1 _ b a s e 6 4 = " y F e 8 s e V d C S C u p 3 8 I r / J y h 1 j A

< l a t e x i t s h a 1 _ b a s e 6 4 = " y F e 8 s e V d C S C u p 3 8 I r / J y h 1 j A

O T T U i 8 1 P K F s T I e 8 a 6 m i E T d + N r 9 0 S s 6 s M i B h r G 0 p J H P 1 9 0 R G I 2 M m U W A 7 I 4 o j s + z N x P + 8 b o r h j Z 8 J l a T I F V s s C l N J M C a z t 8 l A a M 5 Q T i y h T A t 7 K 2 E j q i l D G 0 7 J h u A t v 7 x K W p d V r 1 a 9 u q 9 V 6 r U 8 j i K c w C l c g A f X U I c 7 a E A T G I T w D K / w 5 o y d F + f d + V i 0 F p x 8 5 h j + w P n 8 A T 8 O j S M = < / l a t e x i t > µ 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " + K U m X L i 1 t

y 9 a S V 8 y c w h 9 4 n z + E B I 5 4 < / l a t e x i t > µ 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " B A 9 C X K X A s 9 + O c N 5 8 I O V o 8 c V C W Q A = " > A A A B 7 H i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e y W g n o r e P F Y w W 2 F d i n Z N N u G J t k l y Q p l 6 W / w 4 k E R r / 4 g b / 4 b 0 3 Y P 2 v p g 4 P H e D D P z o l R w Y z 3 v G 5 U 2 N r e 2 d 8 q 7 l b 3 9 g 8 O j 6 v F J x y S Z p i y g

X H e n Y 9 l a 8 H J Z 0 7 h D 5 z P H w K w j Z I = < / l a t e x i t > t 2

< l a t e x i t s h a 1 _ b a s e 6 4 = " q 4 Z 1 m j f / m 0 K

< l a t e x i t s h a 1 _ b a s e 6 4 = " y F e 8 s e V d C S C u p 3 8 I r / J y h 1 j A The correctness of the satisfiability distances introduced above is stated with the following proposition.

Proposition 6.2.1 (Soundness of distances from Definition 6.2.4)

Let σ ∈ P ath(M) be a path of a CTMC M, and φ an MITL temporal that has the form

We should first prove this result for

µ where µ is elementary.

Lemma 6.2.1

For an elementary proposition µ and an MITL temporal formula φ

Estimation of the constant C

An estimation of C can be obtained by a single-point estimation of P r(φ; M θ * ) and π φ-ABC (θ * ). θ * should be chosen wisely: verifying φ should not be rare and θ * should be in a region where π φ-ABC can be well approximated (a region of high posterior probability). Then, P r(φ; M θ * ) can be estimated with statistical model checkers. In order to have a more stable kernel density estimation, one can also choose several θ * and estimate the constants and compute the mean.

Applications

We applied the automaton-ABC method with satisfiability distances to estimate the satisfaction probability function of Poisson processes and three models of biological systems: an enzymatic reaction network (Michaelis-Menten kinetics), a model of viral infection and SIR chemical reaction network. Our results are compared with model checking or statistical model checking from Prism or Cosmos. We also provide an estimation of the experiments' computational times.

An example with Poisson processes

Let M λ be the CTMC defined on a space set S = N by the initial distribution and the transition rate matrix Q: α = δ 0 (the initial state is 0)

We have defined a Poisson process of parameter λ, described by the Chemical Reaction Network

where the kinetic rate is constant and equals λ (it does not follow mass-action law).

The associated stochastic process is denoted S λ = (S λ t ) t≥0 . Let us solve the reachability problem:

Remarks

Results have been obtained by running Algorithm 18 with a sample size N = 1000.

There are no notable differences in performance between Algorithm 17 and 18 for 1D experiments on regions T R1, T R2, T R3 because of the large size of the resulting distributions. However, simple automaton-ABC Algorithm 17 is not worth considering for T R4 and T R5. Given the large size of the support for the considered priors ([0, 100]×[0, 100]) we remark that the probability of sampling (a pair of) parameters in the resulting posterior-distribution is about 90×0.03 2 100 * 100 ≈ 10 -4 . This leads to a tiny probability of drawing from the prior N = 1000 particles that fall in such a narrow distribution, let alone that even a parameter sampled in the obtained distribution could produce paths that do not satisfy φ. By adding several transitional steps with the SMC version (Algorithm 18), the problem becomes treatable. The results for T R4 required about 2 × 10 5 simulations of the model, which is the highest number of simulations in all experiments.

SIR

We consider the SIR compartmental model [START_REF] Kermack | A Contribution to the Mathematical Theory of Epidemics[END_REF] already described in Equations 2.15. We tested our algorithm on a single formula φ = G [0,100] (I > 0) ∧ F [100,120] (I = 0) with one 1D experiment and one 2D experiment. This formula means that the considered epidemic is active within the time interval [0, 100] but disappears during the time interval [100,120]. In the 1D experiment, k i = 0.0012 is fixed, whereas k r ∼ U (0.005, 0.2). In the 2D experiment, both parameters vary: k i ∼ U (5•10 -4 , 0.003) and k r ∼ U (0.005, 0.2). Figure 6.13 reports the results of both experiments, including comparing of the satisfaction probability function obtained with the automaton-ABC method with that obtained with the Prism model checker (numerical method for both 1D and 2D experiments). One can see that the satisfaction probability function is well reconstructed in both cases.

Whereas success is quite expected in the 1D experiment because the histogram suggests a Gaussian-like shape of the density, the result for the 2D experiment is more remarkable because it is hard to guess the density shape based on the 2D histogram.

However, our algorithm managed to reproduce the same complex shape and values given by the Model Checking estimates. [0, 100]. It lasted more than one hour, and we saw this grid missed a region of high probability for φ 4 formula.

Exp

Script 10 (Automaton-ABC experiments related scripts)

Several scripts can be found in the git repository of the thesis:

• code/chap5/prism are related to model checking and statistical model checking estimations of the satisfaction probability function.

• code/chap5/cosmos are related to the computation of average distances with Cosmos.

• code/chap5/julia are related to the run of automaton-ABC algorithm, which are based on our package (Appendix A), as well as illustrations of the Applications Section 6.5.

A comparison with Smoothed Model Checking

In this chapter, we have presented a method that tackles the estimation problem (Definition) in stochastic model checking. As seen in Chapter 3, it is in line with the work of [START_REF] Bortolussi | Smoothed model checking for uncertain Continuous-Time Markov Chains[END_REF] where the so-called Smoothed Model Checking (Smoothed MC) method estimates the satisfaction probability function of parametric CTMC. Methods belonging to this family do not solve the parameter synthesis problem explicitely but provide reasonable estimates of the satisfaction function that can lead to parameter synthesis, e.g. [START_REF] Bortolussi | Bayesian statistical parameter synthesis for linear temporal properties of stochastic models[END_REF] which is also based on Smoothed MC. In this section, we discuss some results of Smoothed MC algorithm for which a Python version is given in [START_REF] Bortolussi | Bayesian statistical parameter synthesis for linear temporal properties of stochastic models[END_REF]. We have already described this method in Section 4.3.1. For the sake of comparison, we 6.8. Summary

Linear Hybrid Automata for non-elementary regions

The presented automata and results were developed for elementary regions µ. However, the automata are extendable to non-elementary regions. Here, we discuss two cases that induce non-elementary regions and describe automatic rules to construct automata for these non-elementary regions.

First, when the condition implies two species O and P , we have to make disjunction elimination. Let us consider µ = (x

2 ) ∨ (x

2 ). We duplicate automaton variables (e.g d becomes d (O) and d (P ) ). Each asynchronous edge of the two automata can be written as

• γ t refers to the conditions relative to the time t and γ O the other ones (these are the conditions that depend on x O , x

• U O refers to the updates of variables related to O (all of the updates are related to O in our case).

Thus, any asynchronous edge has to be replaced by three edges:

Each edge that goes to the final location has to update the final distance, d f inal = min(d (O) , d (P ) ). Also, the synchronised transitions must update n (O) and n (P ) .

Second, for a species O, we consider µ

when it is not an elementary region (it means x 1 ≤ x 2 < x 3 ≤ x 4 ). Thus, one can apply the methodology detailed above by taking P = O and x

These rules can be iterated to obtain automata with non-elementary regions that define n ∈ N elementary regions. Finally, one can construct a new A G∧F automata based on these transformations. Indeed, this automaton is composed of the two automata A G and A F linked by an asynchronous transition.

Summary

In this Chapter: Chapter 6. Automaton-ABC for the statistical parametric verification of CTMCs

145

• We have developed a new method to estimate the satisfaction probability function related to time-bounded reachability problems. The new procedure involves automaton-ABC with satisfiability distances between a trajectory and a property, computed by synchronised simulation of a CTMC with a Linear Hybrid Automaton.

• We have proved their efficiency in terms of correctness of the results and the efficiency of the implemented methods over several biological models.

Several tracks are worth considering to continue and improve this work:

• Any satisfiability distance in accordance with an MITL formula can be considered if the Proposition 6.2.1 still holds.

• Develop a program in our package that automatically creates the extended automata based on rules detailed in Section 6.7.2, and design an automaton for all types of Until formulae.

• Creates a new procedure of parameter synthesis based on automaton-ABC. Indeed, in Algorithm 18 the considered ϵ tolerance is zero in this new ABC procedure: it guarantees the link of the ABC posterior with the satisfaction probability function. However, one can think to relax this tolerance to have better confidence over the subset of parameters that do not fulfil the property φ (i.e. the complementary of the parameter space subset estimated by automaton-ABC).

• Describe how automaton-ABC statistically behaves if the trajectories are approximatively simulated with, for example, Tauleap Algorithm 3. would be interesting for more sophisticated runs of automaton-ABC.

Last words

Many phenomena in Systems Biology involve complex behaviours that are stochastic and complex. Performing inference or verification on such models is generally difficult. In this thesis, we attempted to show how Bayesian inference and statistical verification, via the Approximate Bayesian Computation (ABC) methods and the Hybrid Automata Stochastic Language (HASL), could both benefit from each other. Bayesian inference has benefited from HASL by its capacity of formally describe such behaviours. Statistical verification has benefited from Bayesian inference by its efficiency to explore the parameter space of statistical models.

A.1 Introduction

This thesis introduced a framework for statistical inference and verification of CTMCs with quite heavy notations and concepts. Thus, there is a need of developing an interface that allows users to be familiar with these methods easily and run them efficiently.

We choose Julia [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF] as a programming language to address this task because it is both an efficient language and an easily usable one. Julia is a high-level multi-paradigm language, combining functional and object-oriented programming. It is adapted to mathematical notations, near Python syntax, can easily call and import Python and R objects, and can interact easily with C and Fortran code.

We wanted to construct a package that implements the concepts presented in this thesis with qualities such as:

• Computationally efficient.

• Easy to use.

154

A.2. A few introductory examples

• Deeply tested, both in terms of execution and statistical relevance.

Efforts have been made to conjugate these three requirements. The two first ones, which seem contradictory, are reachable via the use of the Julia language. Indeed, Julia allows a simple mathematical oriented syntax, but performance is reachable with 100% Julia code using computational diagnosis tools mainly described in the web page (JuliaLang, 2020).

Remark A.1.1 (Some Julia notations)

In the following, T1::AT means the type T1 is a subtype of AT . It can be interpreted as "T1 inherits from AT" in classical object-oriented programming. var::T1 means var has type T1 (typeof(var)<:T1 ).

In Julia, naming a function with an exclamation mark (for example, set_x0! ) is a convention to tell the user that one of the variables may be modified to prevent side effects.

Script 11 (MarkovProcesses.jl package)

MarkovProcesses.jl package is available at the git repository https://gitlab-research. centralesupelec.fr/2017bentrioum/markovprocesses.jl.

A.2 A few introductory examples

A.2.1 Simulation of the SIR model

Let us start with the simulation of the SIR model. One can easily create a model that is not available in models/ thanks to the @network_model macro. The parameters and initial state of models created with @network_model are set to zero, and all the variables are observed. One can change that with set_x0! , set_param! , set_time_bound! or set_observed_var! .

A.2.2 Simulation of the ER model synchronised with A F automaton

The package also allows the simulation of a CTMC synchronised with a Linear Hybrid Automaton described in Section 4.4 of Chapter 4. • Line 8 runs the automaton-ABC algorithm and stores the results in res .

One can access the ABC empirical posterior distribution by res.mat_p_end . As there is only one degree of freedom in our parametric model, the size of mat_p_end is 1x1000.

• Line 9 plots the histogram of the ABC empirical distribution. This is the minimal requirement list. It can, and it will be extended to other features. The idea is to present a framework where statistical inference or verification algorithms can be implemented easily with the framework.

A.3. Structure of the package

The source code is organised as follows:

Appendix A. MarkovProcesses.jl : A Julia package for efficient simulation, statistical inference and verification methods of Markov Processes. 161

• core/common.jl: Contains the definition of the common types (struct ) and constructors of the variables of our package. This file regroups all the required type definitions.

• core/model.jl: Contains methods related to models, for example, simulation, access or modification of fields.

• core/trajectory.jl: Contains methods related to trajectories (a trajectory is a type of variable returned by a simulation of a model).

• core/lha.jl: Contains methods related to Linear Hybrid Automata.

• core/network_model.jl: Implementation of @network_model .

• models/: Each file is an implementation of a model. It creates a variable whose type is T <: Model .

• algorithms/: Contains statistical inference and verification algorithms.

• tests/: Contains tests of our package powered by the Test Julia package.

• bench/: Contains benchmarks of several tasks powered by the BenchmarkTools Julia package.

A.4 Type diagram

In this section, we describe the general structure of the project. The package was developed with an object-oriented point of view because it is more suited for the various concepts in this thesis.

A.5 Implementation

A.5.1 Simulation of ContinuousTimeModel

In this section, we describe the simulation of an object model of type ContinuousTimeModel Ensure: A Trajectory object that contains the simulation.

Init: The variables x n+1 , l_t, l_tr are vectors in which the simulation of the next time step is stored. This kind of implementation is used to decrease memory allocation, which gives better performance. That leads to l_t being a 1-length vector for the next time and l tr being a 1-length vector for the next transition/event.

The pseudo-code is simplified. Among other things such as the observation of the model and how the variables are mapped to an index, we have implemented a buffering system for allocating memory at two levels in order to reach better performance.

• The field model.estim_min_states gives an approximative estimation of the number of states. Then, the simulator will allocate a vector of length model.estim_min_states for each variable instead of pushing the values step by step, requiring (approximately) reallocation at each step during execution.

The vectors are resized if there are not enough simulation steps. By default, the value is set to 50.

• The field model.buffer_size is the length of the buffer after reaching model.estim_min_states . Then, the same buffering system described above is executed. The default value is 10.

Even if there is a slight improvement of the performance, the gain is a bit disappointing because by default, Julia allocates and manages arrays ingeniously (for example preallocation).

Appendix A. MarkovProcesses.jl : A Julia package for efficient simulation, statistical inference and verification methods of Markov Processes. 165

A.5.2 Simulation of SynchronizedModel

In this section, we describe the simulation of an object sync_model of type A.6. Use of ABC methods

A.5.3 About trajectories

• times(σ) returns the times of the transitions related to transitions(σ) (for example, the times when a reaction occurred in the SIR model). If a transition equals nothing , then no transition occurred at the time (for any simulation, transitions(σ)[1] == nothing ). If the trajectory is simulated from a bounded model, then the last transition is nothing and vice versa.

An object SynchronizedTrajectory is similar to Trajectory , except it contains the last state of the LHA.

A.5.4 Synchronisation with LHA

The synchronisation of a ContinuousTimeModel with an LHA is mainly embodied by the next_state! method: • l_obs::Vector{T2} is a collection of observations.

• func_dist must have the signature func_dist(l_sim::Vector{T1}, l_obs::Vector{T2})

where T1 is the type of the simulation returned by the model in pm .

The keyword arguments are:

• nbr_particles is the number of particles

• tolerance is the final tolerance ϵ.

• duration_time is a threshold for the time execution. When it is exceeded, the algorithm finishes its current execution and ends.

• bound_sim is a threshold for the number of simulations. When it is exceeded, the algorithm finishes its current execution and ends.

• If save_iterations equals true , the parameters matrix, the associated weights and computed distances are saved at each iteration.

• If init_mat_p (a parameters' matrix), init_weights (a weights' vector) and init_vec_dist (a distances' vector) are specified, the ABC-SMC algorithm begins with these values instead of simulating from the prior.

Examples can be found in Script 4.

A.6.2 ABC with synchronised simulation

Automaton-ABC-SMC Algorithm 16 is run through the method automaton_abc : Several examples with different automata are available:

• A per period automaton in Script 5,

• A ABC,ϵ ABC euclidean distance automaton in Script 6,

• A F , A G , A G∧F automata for satisfiability distances of time-bounded reachability in Script 10.

A.7 Tests

In the following, we describe the tests of the package located in tests/. Each subdirectory represents a specific group of tests, and each Julia file of a sub-directory is a test.

For example, tests/unit contains all the unit tests of the package. Each file tests a specific characteristic. These tests are grouped in the file tests/run_unit.jl with the Test package and @testset macro. Running tests/run_unit.jl runs all the unit tests contained in tests/unit. This is done for every sub-directory.

Appendix A. MarkovProcesses.jl : A Julia package for efficient simulation, statistical inference and verification methods of Markov Processes. 169

One can run tests/run_all.jl to run all the tests group.

A.7.1 Execution test

These tests check if the execution behaves correctly for essential package functions.

• tests/unit groups the unit tests.

• tests/dist_lp groups the tests of the L p distance between two trajectories.

Some trajectories are hard-written to test specific behaviours, whereas other computations are tested with a Riemann sum.

• tests/automata groups tests of linear hybrid automata and synchronised simulations.

• tests/automaton_abc checks if automaton_abc method runs without errors.

• tests/abc_smc checks if abc_smc method runs without errors.

A.7.2 Cosmos based statistical tests

A more original approach to test our methods is done by statistically testing the behaviour of our algorithms with the statistical model checker Cosmos [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF]. Each automaton implemented in our package is written in Cosmos with a .lha file. For now, only the distance values of A F , A G , A G∧F automata are tested.

The basic structure of a test follows:

• A configuration is set:

-Type of model, initial state and parameters.

-Type of automaton and its constants.

-Statistical configuration: width, level

• The test runs a system call of Cosmos and parses the results files.

• The same configuration is run by the package, with the same Cosmos number of simulations.

• It tests if the package result of the average distance is statistically consistent with the width of the confidence interval. It also tests if all the simulated trajectories are accepted by both our package and Cosmos runs because automata A F , A G , A G∧F were designed to accept all the trajectories.

A.8. Benchmarks

These tests are located in tests/cosmos. For now, there are three Julia files, one file per automaton.

• tests/cosmos/distance_F/ER_1D.jl runs the test described above for R1, R2 and R3 configurations of the A F automaton for the ER model described in Section 6.5.2. k 3 varies over the range [0, 100] with a step of 10.

• tests/cosmos/distance_G/ER_R5.jl runs the test for R5 configurations of the automaton A G for ER model. k 1 varies in {0, 0.5, 1.0, 1.5} and k 2 in {0, 40, 80}.

• tests/cosmos/distance_G_F/ER_R6.jl runs the test for the R5 configuration of the automaton A G∧F for the ER model. k 1 varies in {0, 0.5, 1.0, 1.5} and k 2 varies in {0, 40, 80}.

We do not pretend these tests guarantee correct behaviours of any user-defined automaton, but the distances automata can be relatively trusted.

A.8 Benchmarks

In this thesis, the computational efficiency of simulations is critical because likelihoodfree inference and statistical verification methods require a massive amount of simulations. In this section, we present benchmarks of the execution time of simulations.

In this package, efforts have been made to reduce useless computations and limit memory allocation to reach better performances. All the benchmarks were run on a Dell XPS 9370 with CPU Intel i7 8550U @ 1.8GHz x 8 cores.

Script 15 (Benchmarks of MarkovProcesses.jl)

Benchmark scripts are available at code/appendix_pkg/bench. Each file runs a benchmark and generates L A T E X files, which are included in this appendix.

A.8.1 Versus Cosmos

Cosmos is a program developed in C++. It is a well-known compiled language very often used when performance is critical. Julia is a high-level dynamic language, and Julia's compiler produces LLVM native code, a compiler written in C++.

The idea of this benchmark is to show how close the package can be from a program written in C++. For such a scripting language as Julia, it is easy to use syntaxes that allocate too much memory considering our tasks. A simple illustrating Appendix A. MarkovProcesses.jl : A Julia package for efficient simulation, statistical inference and verification methods of Markov Processes. 171 example is Julia's slicing array. The syntax v[m:n] creates a copy of the accessed vector v, which can be too computationally expensive if we only want to read these values (for the record, it is solvable with @view ). Thus, we paid much attention to the computational cost of methods that can be executed a considerable number of times, such as model.f! or next_state! . For example, any ContinuousTimeModel created by @network_model creates a function f! that allocates 0 bytes after a few executions.

We designed several benchmarks:

• Each benchmark implies the following steps:

• A configuration is set:

-Type of model, initial state and parameters.

-Type of automaton and its constants.

-Statistical configuration: width of confidence interval, confidence level interval

• The benchmark runs a system call of Cosmos and parses the results files: number of simulations, time, used memory.

• The same task is run by the package with the same number of simulations. Cosmos. However, the execution times are of the same order (the difference is always less than the execution time itself).

We can observe that the package allocates less memory compared to Cosmos. A deep memory allocation analysis was performed on the package to reduce efficiently useless allocations, which implies some specific knowledge about the Julia language.

These encouraging results need to be put into perspective because our package does not have all the functionalities of Cosmos, which is a more evolved machinery.

Our package does not handle all HASL expressions, which requires computations at each simulation step. For now, only the last value of an automaton variable is available. Also, our benchmarks were not scaled to more complex models. Thus, it is impossible to state that our package will perform better for every model and task.

However, it shows that our package is comparable in terms of execution time and performs better for the tasks considered in the thesis. This is quite challenging for such a high-level dynamic language as Julia compared to a C++ written program.

Appendix Thus, the combination of the two packages can simulate any ContinuousTimeModel , except that there is less flexibility in their @reaction_network macro than ours.

Indeed, in Catalyst.jl, they assume mass-action law for the kinetics, whereas any form can be written for the propensity function in our package. In the following, we can compare the performance of our package to the combination of Catalyst.jl and DifferentialEquations.jl.

We made two benchmarks of performance:

• Benchmark 1: One simulation of the ER model with parameter vector p = [0.2, 40.0, 1.0]. This parameter vector is chosen because the number of steps is significant (at least 7000 states for a simulated trajectory in practice). These benchmarks do not even take into account the loading of modules. For example, a using Catalyst can take about 20 seconds at the first execution of a Julia session (but not the first-ever, which means the package is already precompiled), whereas using MarkovProcesses takes about 3 seconds. This is non-negligible if one distributes computations over an HPC cluster, for example.

A.9 Conclusion

A.9.1 Summary

We have developed a Julia package to simulate CTMCs and synchronised CTMCs with a Linear Hybrid Automata, and the automaton-based ABC methods developed in this thesis. The call of methods is high-level (it takes only a few lines). Tests were developed to verify our experiments statistically. Benchmarks were made to compare performance to other existing programs. We showed that our package performs better than one of Julia's reference package. Execution times are of the same order, which is encouraging, considering that Cosmos is a C++ program.

A.9.2 Perspectives

However, we believe this work builds the foundations of a tool with room for improvement in many ways.

• Use of any form of HASL trajectory expression in automaton-ABC.

• A macro that facilitates the creation of LHA. For now, they are Julia handwritten coded.

Appendix A. MarkovProcesses.jl : A Julia package for efficient simulation, statistical inference and verification methods of Markov Processes. 175

• Tests with other more complicated models.

• Statistical tests of created models with --loop option of Cosmos.

• Statistical tests with automata that are not related to time-bounded reachability.

For the integration of this package into Julia's community, perspectives are:

• Creation of a ReactionSystem object of the ModelingToolkit.jl API.

• Analyse how our methods can be integrated into the Turing.jl framework, which regroups many Bayesian inference methods.

Appendix B

Measure theory

This appendix aims at giving the main tools of probability theory for the construction of probability measures.

Section 2.2 constructs a measure and proves the equality of measures over a σ-algebra generated by a semiring. The consistency of the construction of these measures is based on Caratheodory's extension theorem.

B.1 Results from measure theory: Caratheodory's theorem

The definitions and theorems presented in this section are based on [START_REF] Klenke | Probability Theory[END_REF] Chapter 1).

Let X be a set of elements. In the following definitions, A ⊂ P(X) is a subset of the power set of X.

Definition B.1.1 (Semiring of sets.)

Let X be a set of elements and A ⊂ P(X). A is called a semiring of sets if :

Let µ : R → R ≥0 be a content (an additive set function) over a semiring R. µ is

The outer measure µ * : P(X) → R ≥0 of a set function µ over a semiring R is defined as :

It is monotone and σ-subadditive.

Theorem B.1.1 (Caratheodory's theorem for extended measure)

Let R be a semiring and let µ 0 : R → R ≥0 . If µ 0 is :

Then it admits a unique extended measure µ defined on F(R), and µ 0 = µ |R = µ * |R .

B.2 Semiring of sets of P ath(M)

Let M a CTMC (Section 2.1.2). The probability measure P r M is defined on the cylinder sets: 179

Let R be the set that contains all the cylinder sets C(s 0 , I 0 , . . . , I k-1 , s k ) (R ⊂ P(P ath(M)).

Proposition B.2.1 (Semiring of sets of trajectories.)

Let M be a CTMC defined by an initial distribution α and a transition rate matrix

and for i ∈ {0, . . . , k -1},

Otherwise, let us distinguish two cases.

In other words, A \ B is the set of paths in A that does not have the sequence s This union is disjoint and also finite if for each state the number of transitions is finite.

iii) Stable by

In both cases A ∩ B ∈ R.

B.3 Statistical model of CTMCs

B.3.1 Density of a CTMC

In classical statistical theory, we often consider a set of probability measures (P θ ) θ∈Θ called a statistical model. To have good statistical guarantees, we want our family of measures to be dominated by some measure. In this section, we construct a measure µ independent of P S so that µ dominates P S (i.e P S ≪ µ, absolute continuity).

The set of trajectories of a CTMC can be seen as a subset of (S, R ≥0 ) N . Then, C(k, s 0:k , I 0:k-1 ) can be rewritten as {(s 0 , . . . , s k )} × (I 0:k-1 ∪ {+∞}) which is a subset of (S, R ≥0 ) k+1 .

181

We can define a product measure on k i=0 (S, R ≥0 ). δ k+1 is the counting measure on S k+1 and λ k is the Lebesgue measure on R k .

Let k ∈ N * . Let µ : R → R ≥0 be the monotone set function defined on the semiring R :

Basically, this operation means that if A ∈ R, then the measure of A is the volume of I 0 × . . . × I k-1 . The measure of any intersection of sets of R is defined because it is a semiring, so stable by ∩.

The measure of any union of sets in F(R) is defined by the outer measure µ * in Definition B.1.5.

Remark B.3.1

On any set A = C(s 0 , I 0 , . . . , I k-1 , s k ), A can be represented as {s 0:k } × I 0 × . . . × I k-1 .

In this case, µ(A) = (δ (k+1) ⊗ λ Proof. (i) (ii) Here we do not dive into the details, but the σ-subadditivity of µ * becomes an additivity in F(R). This is done by the definition of M(µ) that defines the µ-measurable sets which is proven to be a σ-algebra, so F(R) is contained in it (Lemma 1.50 and 1.51 in [START_REF] Klenke | Probability Theory[END_REF]. i.e. a stationary sequence of integers. The set of stationary sequences of integers called S 0 is countable (as a countable union of the sets of stationary sequences of order i ∈ N).

{C(s 0:k , [0, n] k )} (s 0:k ,n)∈S 0 is countable, the union of all is (S, R ≥0 ) N (set of any trajectory) and ∀(s 0:k , n) ∈ S 0 , µ(C(s 0:k , [0, n])) < +∞ (by definition). µ is σ-finite. We have a series with positive values whose limit sum is equal to zero, hence all the terms of the series are equal to zero.

On the other hand, µ(A n ) = λ kn (I (An) 0:kn-1 ) = 0, i.e. ∀n, ∃i < k n , ]a

i=0 (e -E(s i )a (An)

i -e -E(s i )b (An) i ) = 0.

Then P S * (A) ≤ +∞ n=0 P S (A n ) = 0 by definition of the outer measure with inf.

We can conclude P S ≪ µ.

Corollary B.3.1 (Density of a CTMC)

Let S be a CTMC. P S admits a density p S with respect to µ :

Proof. Let S be a time-homogeneous CTMC with a infinitesimal generator matrix Then, ∀σ = s 0 t 0 -→ . . .

---→ s k , p S (σ) = P (s 0 ) k i=1 Q(s i-1 , s i )e -E(s i-1 )t i-1 .

Appendix C

A simple analytical example of

ABC inference

Let us illustrate the ABC method and kernel density estimation with a simple example. We consider n i.i.d observations y exp = (y (1) , . . . , y (n) ) from N (µ, σ 2 0 ). We want to estimate θ = µ via Approximate Bayesian Computation (σ 0 = 1.0 is known). The prior is uniform over [-1, 1] (p(θ) ∝ 1).

C.1 Computation of the true posterior and ABC posterior

In this case, the different posteriors can be computed analytically, and the mean estimator is a sufficient statistic. We denote η(y exp ) = y exp = 1 n n i=1 y (i) .

First, let us compute the true posterior distribution.

C.2 Simulations

We run the ABC-SMC algorithm with 100, 1000 and 10000 particles.

With ϵ = 0.01. Les trois premiers chapitres réfèrent à l'état de l'art des différents domaines