Alexandre David

Kim Guldstrand Larsen

Axel Legay

Marius Mikučionis

Danny Bøgsted Poulsen

Sean Sedwards

Runtime Verification of Biological Systems

Keywords: runtime verification, synthetic biology, statistical model checking, genetic oscillator, MITL, frequency domain analysis, Uppaalsmc, Plasma

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

Introduction

It is conceivable to design systems in such a way that makes their analysis easier, but it is most usually the case that they are optimised for other constraints (efficiency, size, cost, etc.) and that they evolve over time, developing highly complex and unforeseen interactions and redundancies. These phenomena are epitomised by biological systems, which have absolutely no inherent need to be understandable or analysable. The discovery that the genetic recipe of life is written with just four characters (nucleotides Adenine, Cytosine, Guanine and Thymine) that are algorithmically transcribed and translated into the machinery of the cell (RNA and proteins) has led scientists to believe that biology also works in a computational way. The further realisation that biological molecules and interactions are discrete and stochastic then suggests the idea that biological systems may be analysed using the same tools used to verify, say, a complex aircraft control system.

Using formal methods to investigate natural systems can thus be seen as a way to challenge and refine the process of investigating man-made systems. It is very difficult to reason about systems of this type at the level of their descriptions, however. It is much more convenient to directly analyse their observed behaviour. In the context of computational systems we refer to this approach as runtime verification, while in the case of biological systems this generally takes the form of monitoring the simulation traces of executable computational models.

To demonstrate runtime verification of biological systems we apply advanced statistical model checking (SMC) techniques to a model of robust biological oscillations taken from the literature. SMC works by verifying multiple independent simulation traces of a probabilistic model against a property specified in linear temporal logic. The results are then used in an hypothesis test or to estimate the probability of the property. In adopting this approach, SMC avoids constructing the generally intractable explicit representation of the state space of the system. The price paid is that results are only known within confidence intervals, however these may be made arbitrarily tight by increasing the number of simulation runs. SMC can thus be seen as a specific instance of runtime verification.

The model we have chosen to investigate demonstrates some of the mechanisms used by biological systems to maintain reliable performance in the face of inherent stochasticity. These mechanisms are literally vital and have relevance beyond biology (e.g. amorphous computing). To perform our investigation we use two recently developed SMC platforms: that incorporated in Uppaal and Plasma. Uppaal-smc offers a generic modelling language based on stochastic hybrid automata, while Plasma aims at domain specific support and here accepts biological models represented in chemical syntax. Although our chosen model was conceived to be stochastic, its original description and analysis were in the continuous (ODE) domain. We therefore compare the behaviour of deterministic and stochastic models by performing a frequency domain analysis, taking advantage of Uppaal's recently implemented ability to work with ODE representations. We verify various interesting temporal properties of the model and compare Plasma's direct implementation of bounded LTL with Uppaal's monitor-and rewrite-based implementations of weighted MITL.

Beyond Runtime Verification with SMC

Runtime verification (RV) [START_REF] Havelund | Synthesizing monitors for safety properties[END_REF][START_REF] Rosu | Allen linear (interval) temporal logic -translation to ltl and monitor synthesis[END_REF] refers to a series of techniques whose main objective is to instrument the specification of a system (code, ...) in order to disprove potentially complex properties at the execution level. The main problem of the runtime verification approach is that it does not permit to assess the overall correctness of the entire system.

Statistical model checking (SMC) [START_REF] Bensalem | Statistical model checking: Present and future[END_REF][START_REF] Younes | YMER: A Statistical Model Checker[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF] extends runtime verification capabilities by exploiting statistical algorithms in order to get some evidence that a given system satisfies some property. Given a program B and a trace-based property3 φ , Statistical model checking refers to a series of simulation-based techniques that can be used to answer two questions: (1) Qualitative: is the probability for B to satisfy φ greater or equal to a certain threshold θ (or greater or equal to the probability to satisfy another property φ ′ )? and (2) Quantitative: what is the probability for B to satisfy φ?

We briefly review SMC approaches, referring the reader to [START_REF] Bensalem | Statistical model checking: Present and future[END_REF][START_REF] Younes | YMER: A Statistical Model Checker[END_REF] for more details. The main approaches [START_REF] Younes | YMER: A Statistical Model Checker[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF] proposed to answer the qualitative question are based on hypothesis testing. Let p be the probability of B |= φ, to determine whether p ≥ θ, we can test H : p ≥ θ against K : p < θ. A test-based solution does not guarantee a correct result but it is possible to bound the probability of making an error. The strength (α, β) of a test is determined by two parameters, α and β, such that the probability of accepting K (respectively, H) when H (respectively, K) holds is less or equal to α (respectively, β). Since it impossible to ensure a low probability for both types of errors simultaneously (see [START_REF] Younes | YMER: A Statistical Model Checker[END_REF] for details), a solution is to use an indifference region [p 1 , p 0 ] (with θ in [p 1 , p 0 ]) and to test H 0 : p ≥ p 0 against H 1 : p ≤ p 1 . Several hypothesis testing algorithms exist in the literature. Younes [START_REF] Younes | YMER: A Statistical Model Checker[END_REF] proposed a logarithmic based algorithm that given p 0 , p 1 , α and β implements the Sequential Ratio Testing Procedure (SPRT) (see [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF] for details). When one has to test θ≥1 or θ≥0, it is however better to use Single Sampling Plan (SSP) (see [START_REF] Younes | YMER: A Statistical Model Checker[END_REF][START_REF] Bensalem | Statistical model checking: Present and future[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF] for details) that is another algorithm whose number of simulations is pre-computed in advance. In general, this number is higher than the one needed by SPRT, but is known to be optimal for the above mentioned values. More details about hypothesis testing algorithms and a comparison between SSP and SPRT can be found in [START_REF] Bensalem | Statistical model checking: Present and future[END_REF].

In [START_REF] Hérault | Approximate probabilistic model checking[END_REF][START_REF] Laplante | Probabilistic abstraction for model checking: An approach based on property testing[END_REF] Peyronnet et al. propose an estimation procedure (PESTIMA-TION) to compute the probability p for B to satisfy φ. Given a precision δ, Peyronnet's procedure computes a value for p ′ such that |p ′ -p|≤δ with confidence 1 -α. The procedure is based on the Chernoff-Hoeffding bound [START_REF] Hoeffding | Probability inequalities[END_REF].

Model and Properties

A genetic oscillator

It is well accepted that molecules are discrete entities and that molecular interactions are discrete. It is further accepted that molecules move randomly as a result of collisions with other molecules (thermal noise). From this it can be inferred that chemical reactions are the result of random interactions and can be modelled as stochastic processes. Biological organisms based on chemical reactions are thus supreme examples of complex stochastic systems. The means by which they overcome low level non-determinism and achieve apparent high level determinism is the subject of much ongoing research and informs such fields as amorphous computing [START_REF] Abelson | Amorphous computing[END_REF].

Oscillation, arising from the execution loops in computer programs, is of great relevance to runtime verification of automated systems. Oscillation also plays a crucial role in biology -life being an essentially cyclic process. One of the key oscillatory behaviours in biology is the circadian rhythm that allows an organism to take advantage of periods of day and night to optimise when to maximise activity and recovery. In light of this, we have chosen the genetic circadian oscillator of [START_REF] Barkai | Biological rhythms: Circadian clocks limited by noise[END_REF][START_REF] Vilar | Mechanisms of noiseresistance in genetic oscillators[END_REF] as the focus of our analysis. This synthetic model distils the essence of several real circadian oscillators and demonstrates how a reliable system can be constructed in the face of inherent stochasticity. In particular, the model has been shown in [START_REF] Hilborn | Stochastic coherence in an oscillatory gene circuit model[END_REF] to exhibit a kind of regularity referred to as stochastic coherence.

Though the authors of the model were interested in its stochastic properties, they nevertheless chose to represent it in the form of a system of ordinary differential equations (ODEs, reproduced in Figure 1). Each of the equations describes the infinitesimal rate of change of a particular molecular species; the functions being the sums of the rates of all reactions involving the species, weighted by the direction (positive or negative, corresponding to creation and consumption) and size of the corresponding change. ODEs are commonly used to represent the dynamics of chemically reacting systems and it is traditional to consider concentrations (numbers of molecules per unit volume). Trajectories of ODEs can closely approximate stochastic dynamics when the system operates near to the thermodynamic limit (infinite population sizes). This is rarely the case with biological models of cellular processes, which frequently consider molecular species in very low copy numbers. An obvious example is that within a cell there is often just a single copy of a particular gene (as in the case of the genetic oscillator we describe here). The ODE trajectory is often considered (informally) to be the 'average' of the stochastic traces, implying that the noise is somehow superimposed on top of a deterministic trajectory. In fact, the noise is an inherent part of the stochastic trajectory and the ODE describes the behaviour of the limit of the stochastic process as populations are taken to infinity while maintaining the same concentrations [START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF]. We demonstrate this using frequency domain analysis in Section 3.4.

Using a standard translation between deterministic and stochastic semantics of chemically reacting systems (see, e.g., [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]), it is possible to transform the ODEs given in Figure 1 into the chemical reaction syntax of Equations [START_REF] Abelson | Amorphous computing[END_REF][START_REF] Alur | A theory of timed automata[END_REF][START_REF] Barkai | Biological rhythms: Circadian clocks limited by noise[END_REF][START_REF] Bensalem | Statistical model checking: Present and future[END_REF][START_REF] Bulychev | Monitor-based statistical model checking for weighted metric temporal logic[END_REF][START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF][START_REF] Havelund | Synthesizing monitors for safety properties[END_REF][START_REF] Hérault | Approximate probabilistic model checking[END_REF][START_REF] Hilborn | Stochastic coherence in an oscillatory gene circuit model[END_REF][START_REF] Hoeffding | Probability inequalities[END_REF][START_REF] Ihekwaba | Communicating oscillatory networks: frequency domain analysis[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF]. These can then be visualised as a stochastic Petri net (Figure 2). The model comprises two genes (D A and D R ) that are transcribed (Equations (5-8)) to produce two micro-RNA molecules (M A and M R , respectively) that are translated (Equations (9,10)) to produce two proteins (A and R, respectively). A acts as a promoter for its own gene (Equation ( 1)) and for that of R (Equation (3)) by reacting with D A and D R to produce their more efficient active forms D ′ A and D ′ R . A and R dimerise (Equation ( 11)) to form complex protein C that eventually degrades back to D. Oscillation arises from the fact that A is part of a positive feedback loop involved in its own production and promotes the production of R that, in turn, sequesters A (i.e., removes it) via the production of C (Equation ( 11)). This mechanism and other mechanisms of biological oscillation are discussed in more detail in a recent review of synthetic oscillators [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF].

dDA/dt = θAD ′ A -γADAA dDR/dt = θRD ′ R -γRDRA dD ′ A /dt = γADAA -θAD ′ A dD ′ R /dt = γRDRA -θRD ′ R dMA/dt = α ′ A D ′ A + αADA -δM A MA dMR/dt = α ′ R D ′ R + αRDR -δM R MR dA/dt = βAMA + θAD ′ A + θRD ′ R -A(γADA + γRDR + γC R + δA) dR/dt = βRMR -γC AR + δAC -δRR dC/dt = γC AR -δAC
A + DA γ A -→ D ′ A (1)
D ′ A θ A -→ DA + A (2) A + DR α R -→ D ′ R (3)
D ′ R θ R -→ DR + A (4) D ′ A α ′ A -→ MA + D ′ A (5) DA α A -→ MA + DA (6) D ′ R α ′ R -→ MR + D ′ R (7) DR α R -→ MR + DR (8) MA β A -→ MA + A (9) MR β R -→ MR + R (10) A + R γ C -→ C (11) C δ A -→ R (12) A δ A -→ ∅ (13) R δ R -→ ∅ (14)
MA

δ M A -→ ∅ (15)
MR

δ M R -→ ∅ (16)
Each equation describes a possible productive interaction between types of molecules (molecular species). Monomolecular reactions of the form A k → • • • have the semantics that a molecule of type A will spontaneously decay to some product(s) following a time delay drawn from an exponential distribution with mean k. Bimolecular reactions of the form A + B k → • • • have the semantics that if a molecule of type A encounters a molecule of type B they will react to become some product(s) following a time delay drawn from an exponential distribution with mean k. It is usually assumed that the system is 'well stirred' [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF], such that the probability of a molecule occupying a particular position is uniform over the physical space of the system. This rarely represents reality in the case of biological cells, however it is a widely used mathematical expedient that is common to both deterministic and stochastic treatments. The consequence is that molecules lose their individuality (populations are treated as multisets) and the rate of molecular interactions is simply the product of a rate constant (k in the examples) and a combinatorial factor arising from population sizes [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF]. This is known as mass action kinetics [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF]. Referring to molecular species A and B in the example reactions given above, for populations of instantaneous size A and B the overall rates of monomolecular and bimolecular reactions are given by kA and kAB, respectively.

Properties

The language we use to describe properties is based on the dense timed logic MITL, having the form

φ = φ ∨ φ | φ ∧ φ | ¬φ | ⊤ | ⊥ | φU [a;b] φ | φR [a;b] φ | Xφ | α

Properties of the Oscillator

The model exhibit an oscillatoric behaviour in which a peak of one protein is followed by the increase of another protein. The increase of one protein also appears to be governed by highly regularity in the sense that one peak level is followed by another peak level in a specific amount of time.

In order to detect peaks we first define the shape of a peak. We say there is a peak if the protein level is above a threshold thres L and within l time units drops below another threshold thres R .

Using M IT L ≤ we can express that we at the given time is in a peak of the N variable as

φ peakN ≡ N > thres L ∧ ♦ [0;l] N < thres R .
Expressing that there is a periodicity in the peaks of a single variable N within the first 1000 time units can be done using the formula:

≤1000 (φ peakN =⇒ ♦ ≤p φ peakN ),
where p is the maximum time between peaks. The same form of expression can of course also be used to express that a peak on the N variable should be followed by a peak on the M variable.

Frequency domain analysis

Frequency domain analysis provides a rigorous yet intuitive means to quantify the behaviour of stochastic systems from observations of their executions. This methodology is particularly relevant for oscillatory biological systems [START_REF] Ihekwaba | Communicating oscillatory networks: frequency domain analysis[END_REF], but is not limited to these and is able to characterise the distance in behaviour between different models, different systems and different parts within the same system. It can also measure the difference between different simulation algorithms or semantics applied to the same system.

Our technique is to generate N simulation traces sampled at constant time intervals, δt, resulting in K sampled points per simulation. From each set of sampled points and for each state variable of interest we calculate a complex frequency spectrum using a 'fast Fourier transform' (FFT) algorithm. From these we generate N corresponding magnitude spectra and then calculate the pointwise average magnitude spectrum. The average magnitude spectrum often gives a visually compact notion of the complex stochastic behaviour of the system and can also be used to quantify a distance between behaviours using standard statistical metrics.

K and δt are chosen according to the temporal characteristics of the phenomenon of interest: Kδt is the maximum observed simulated time; (Kδt) -1 is the low frequency resolution (the spacing between spectral components) and (2δt) -1 is the maximum observable frequency. It is generally desirable to increase K and reduce δt, but note that an optimal value of δt is usually significantly greater than the minimum time between successive update events, since these often do not apply to the same variable and the highest part of the spectrum is often uninformative. N is chosen according to the stochasticity of the system in relation to the desired discrimination of the metric; large N being desirable.

UPPAAL-SMC

The verification tool Uppaal [START_REF] Larsen | Uppaal in a nutshell[END_REF] provide support for modeling and efficient analysis of real-time systems modeled as networks of timed automata [START_REF] Alur | A theory of timed automata[END_REF]. To ease modeling, the tool comes equiped with a user-friendly GUI for defining and simulating models. Also, the modelling formalism extends basic timed automata with discrete variable over basic, structured and user-defined types that may be modified by user-defined functions written in a Uppaal specific C-like imperative language. The specification language of Uppaal is a fragment of TCTL supporting a range of safety, liveness and bounded liveness properties.

Uppaal-smc is a recent branch of Uppaal which support statistical model checking of stochastic hybrid systems, based on a natural stochastic semantics. Uppaal-smc extends the basic timed automata formalism of Uppaal by allowing rates of clocks to be defined by general expressions possibly including clocks, thus effectively defining ODEs. An overview of the architecture is given in Figure 3. The GUI of the tool allows the user to draw automata templates in the editor, instantiate and compose these into a system, simulate the system for easy validation, verify queries, and visualize quantitative answers in the form of plots in the plot composer. The execution engine of Uppaal-smc implements the stochast semantics of interacting hybrid automata, and includes a proprietary virtual machine for the execution of imperative code of the model.

The specification formalism of Uppaal-smc is that of (weighted) MITL, with respect to which four different statistical model checking components are offered: hypothesis testing, probability estimation, probability comparison and simulation. Here the user may control the accuracy of the analysis by a number of statistical parameters (size of confidence interval, significance level, etc.). Uppaal-smc also provides distributed implementations of the hypothesis testing and probability estimation demonstrating linear speed-up [START_REF] Bulychev | Monitor-based statistical model checking for weighted metric temporal logic[END_REF].

The results generated by the analyses can be processed for visualization in various ways: Gannt charts for individual runs monitoring desired variables, plots of density functions and accumulated distribution functions with respect to given (W)MITL properties. Typically the simulation results in gigabytes of data which are filtered on-the-fly to plot only the relevant points. Execution Engine Server Fig. 3: Architecture of Uppaal-smc.

Graphical

Modeling and Checking in UPPAAL-SMC

A Bouncing Ball Example. To illustrate the expressive power of the stochastic hybrid automata language supported by Uppaal-smc, we consider a simple, yet interesting variant of a bouncing ball. Figure 4(a) gives the principle of a ball bouncing on a floor and being hit by a piston. The hybrid model of the ball is given in Fig. 4(b) where three cases are visible: (i) it can be hit while going up (v<=0), (ii) hit while going down (v<0), or (iii) it bounces on the floor. The invariant on the location describes the trajectory of the ball in the form of two differential equations (v' and p'). The piston in Fig. 4(c) can hit the ball only if its position is high enough (p>=6). The ball will rebound with a random dampening coefficient both on the floor and the piston (given by the random function). The delays between hits of the piston are chosen stochastically according to an exponential distribution of rate 5/2. Semantically, the effect of ODE expressions is achieved by an implicit auxiliary process integrating the values based on a given fixed time step and thus directly competing with the rest of processes when the rest of invariant and guard expressions are evaluated. Five different trajectories of the ball are obtained by the query simulate 5 [<=20]{p} and shown in Fig. 4(d). We may also ask the model-checker to estimate the probability that the ball is still bouncing above the height of 4 after 12 time units with the query:

Pr[<=20](<>(time >= 12 && p >= 4))

Here <>(time>=12 && p>=4) is the Uppaal-smc syntax for the MITL property ψ = ♦(time ≤ 12 ∧ p ≥ 4) and Pr[<=20]ψ denote the probability π that ψ will hold within 20 time-units for a random run. Given this query, the interval [0.152, 0.163] is returned as an estimate for π with confidence 99.9% after generating 152020 runs. We can also test the hypothesis: Pr[<=20](<> time >= 12 && p >= 4) >= 0.15 with a region of indifference of ±0.005 and level of significance of 0.1% after generating 18543 runs.

For the analysis of more general MITL properties properties, Uppaal-smc generates monitoring automata to be put in parallel with the system. Statistical model checking requires that these monitors are determinstistic timed automata. Unfortunately, not all MITL properties may be monitored by determinstic timed automata. Thus Uppaal-smc offers a safe confidence interval based on two monitors corresponding to under-and over-approximations of the set of runs satisfying the particular formula [START_REF] Bulychev | Monitor-based statistical model checking for weighted metric temporal logic[END_REF]. Experimental results have shown that we obtain an exact monitor, and most recently this method has been replace by an exact rewrite technique. Figure 5a shows a snippet from an ODE model of the genetic oscilator, where the coefficients (gammaR, betaA, deltaA) and variables (A, C, R, DR) are initialized with the first urgent transition and then the trajectories are computed based on the ODEs (the last three equations from Fig. 1). A snippet from stochastic genetic oscilator model is shown in Fig. 5b, where each reaction (from Eq. 1, 9 and 12) is modeled by a separate automaton. For example, the first automaton can be read as follows: reaction requires positive quantities of A and DA (guard conditions), one of each is consumed (A--and DA--), and one D A is produced (D A++) with an exponential rate gammaA times the available quantities of A and DA.

PLASMA

Plasma is designed to be a high performance and flexible platform for statistical model checking, able to work with multiple modelling and property specification languages. Its basic architecture is shown in Figure 6 and comprises a user interface, a simulation management module, a virtual machine and modules that compile the model and property specifications. Models and properties are

Modelling Language Module

Specification Language Module

Fig. 6: The architecture of Plasma compiled into proprietary bytecode and executed on Plasma's inbuilt virtual machine. Overall control of the verification process is maintained by the simulation management module according to the options specified by the user. The simulation management module contains various statistical model checking algorithms that implement confidence bounds, such as the Chernoff bound of [START_REF] Hérault | Approximate probabilistic model checking[END_REF] and the sequential hypothesis test of [START_REF] Younes | Probabilistic verification of discrete event systems using acceptance sampling[END_REF], plus an importance sampling engine [START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF] that improves the performance when simulating rare properties. For simulating discrete and continuous time Markov models the virtual machine uses the 'method of arbitrary partial propensities' (MAPP [START_REF] Sedwards | A Natural Computation Approach To Biology: Modelling Cellular Processes and Populations of Cells With Stochastic Models of P Systems[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF]) that is an optimised version of the 'direct method' of [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]. The simulation management module executes the property bytecode that, in turn, executes the model bytecode until sufficient simulation steps are generated. In this way simulation traces contain the minimum number of states necessary to decide a property and the simulation management ensures that the minimum number of simulation traces are requested of the simulator.

Modeling and Checking in Plasma

Modelling languages are built on an underlying semantics of guarded commands [START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF] extended with stochastic rates to resolve non-determinism. These have the form (guard, rate, action), where guard is a logical predicate over the state of the system which enables the command, action updates the state of the system and rate is a function over the state of the system that returns the stochastic rate at which an enabled command is executed. This semantics is equally applicable when the rate is actually a probability and time plays no part.

Plasma is designed to be language neutral, so for the present investigation Plasma adopts a simple chemical syntax modelling language that closely mirrors the style of Equations [START_REF] Abelson | Amorphous computing[END_REF][START_REF] Alur | A theory of timed automata[END_REF][START_REF] Barkai | Biological rhythms: Circadian clocks limited by noise[END_REF][START_REF] Bensalem | Statistical model checking: Present and future[END_REF][START_REF] Bulychev | Monitor-based statistical model checking for weighted metric temporal logic[END_REF][START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF][START_REF] Havelund | Synthesizing monitors for safety properties[END_REF][START_REF] Hérault | Approximate probabilistic model checking[END_REF][START_REF] Hilborn | Stochastic coherence in an oscillatory gene circuit model[END_REF][START_REF] Hoeffding | Probability inequalities[END_REF][START_REF] Ihekwaba | Communicating oscillatory networks: frequency domain analysis[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF]. The structure of the model file follows the form: constant initialisations, species initialisations, list of reactions. In the present context there is an implicit assumption of mass action kinetics [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF][START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF][START_REF] Gillespie | Deterministic limit of stochastic chemical kinetics[END_REF] and rate specifies the mean of an exponential random variable that models the time between successive reaction events; non-determinism being resolved by races between realisations of the random variables of competing reactions. Reactions of the abstract form A + B Plasma verifies properties specified in bounded linear temporal logic of the kind described in Section 3.2. The logic accepts arbitrarily nested temporal formulae and Plasma achieves this using a buffer to store sequences of values of the variables of interest. When formulae are not nested, no buffer is required. Algorithms 1 and 2 illustrate the basic notions of checking non-nested temporal formulae, employing discrete time for clarity. Algorithm 3 is a naive implementation of a nested formula to illustrate the purpose of the buffer and how we improve efficiency. Algorithms 1 and 2 generate and consider states in turn, returning a result as soon as φ is satisfied or not satisfied, respectively. These algorithms store nothing and generate the minimum number of states necessary. Algorithm 3 also only generates new states as required, but since the inner loop requires states further into the future than the outer loop, states are stored by the inner loop for subsequent use by the outer loop. As written, Algorithm 3 is naive because the inner loop re-checks states that it has checked on previous iterations of the outer loop. Plasma therefore records where the decision on the previous iteration was made and then needs only check the states after that. The case with continuous time is more complex because the length of the buffer is not known a priori (there may be an arbitrary number of steps to achieve a given time bound). Plasma overcomes this by creating an initial buffer and then extends it as required.

Rare Events

Rare events pose a challenge to simulation-based approaches, so Plasma includes an importance sampling engine that makes it possible to estimate the probability of a rare property by simulating under a distribution that makes the property less rare. Given a property φ, with true probability γ under distribution P , the standard Monte Carlo estimator of γ is given by γ = 1 N N i=1 z(ω i ), where ω i is the trace of a simulation made underP and z(ω) ∈ {0, 1} indicates whether ω |= φ. In general, N must be chosen significantly greater than 1 γ to accurately estimate γ, hence this is computationally expensive when γ is small. By contrast, the importance sampling estimator is given by γ = 1 N N i=1 z(ω i ) P (ωi) Q(ωi) , where Q is ideally a distribution under which traces that satisfy φ are uniformly more likely and ω i is now the trace of a simulation performed under Q. P Q is called the likelihood ratio and in a discrete event simulation can usually be calculated on the fly in constant time. Since Q is chosen to reproduce φ more frequently, N may be significantly less than 1 γ . The effectiveness of importance sampling relies on finding a suitable Q.

An optimal importance sampling distribution is one under which traces that satisfy the rare property are uniformly more likely, to the exclusion of all traces that do not satisfy the property. It is possible to find such distributions by individually modifying all the transition probabilities in the system [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF], however this is often intractable. Plasma thus parametrises the distribution with a low dimensional vector of parameters applied to the rates of its guarded commands [START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF]. In the case of biological systems of the type considered here, this parametrisation corresponds to the rate constants of reactions.

To demonstrate the application of importance sampling to biological systems we consider a simple chemical system comprising A + B . By multiplying the rate constants k 1 , k 2 , k 3 by importance sampling parameters λ 1 = 1.16, λ 2 = 1.15, λ 3 = 0.69, respectively, Plasma is able to estimate this probability using only 1000 simulation runs (of these approximately 600 satisfy the property). The parameters were generated using the cross-entropy algorithm described in [START_REF] Jegourel | Cross-entropy optimisation of importance sampling parameters for statistical model checking[END_REF].

Experiments

Our first level of validation is to inspect the simulation traces to verify that they are sensible. The result from Uppaal-smc is displayed in Figure 7 where the ODE model yields the same (deterministic) trajectories of a pattern which repeats every 26.2 hours. The stochastic model yields an apparently similar but "noisy" pattern where the amplitude and periodicity are also varying. This intuitive similarity is made more formal by the frequency analysis in Section 6.1. Notice that the signal C starts at zero which allows A to reach higher amplitude ated by a MITL formula: true U[<=1000] (A>1100 & true U[<=5] A<=1000), which is false unless the signal A rises above 1100 and then falls below 1000 within 5 time units. By coordinating two instances of such monitors and letting a stopwatch x run only between the two peaks we estimate the distribution of x as a period estimate. For signal C (R) we register a peak when the signal falls from 2000 to 1900 (1500 to 1400 resp.). The query then is a simple probability estimate whether the second peak is reached. Figure 9 shows estimated period with corresponding values at probability peaks. The most instances are situated around the peak value, but there are other tiny bumps on sides. The bump near zero corresponds to false positive recognition of a peak when the signal includes a saw tooth (common in stochastic simulations). The bump around 48 hours corresponds to a missed peak at 24, due to the local amplitude being too low. Interestingly there is another tiny bump at around 73 in Fig. 9b which correspond to two missed peaks in a row.

A period log prob. density Similarly we can estimate the phase difference between the signals as a distance between peaks of different signals.

Figure 10 shows a probability density distribution of a phase difference between A and R signals, which implies that a peak in A typically leads to a peak in R within 6.21, but there might be one missed peak.

Frequency domain analysis

Figures 11 and 12 compare the average frequency spectra of two variables from our deterministic and stochastic simulations. We observe that the deterministic spectra (black) tend to have sharp, well defined, peaks with discernible harmonics at high frequencies. In contrast, the spectra of stochastic simulations (red) tend to have softened peaks, few discernible high harmonics and contain an apparent continuum of frequencies. This is reflected in the reduced amplitude of these spectra relative to their deterministic counterparts: the amplitude of the original trace is effectively divided amongst many more frequencies. We note that the first three harmonics of the deterministic and stochastic spectra appear to coincide, confirming our expectation and intuition that the the two behaviours are 'similar'.

The apparent thickness of the red lines reflects the fact that the lines describe average spectra generated from stochastic data. Increasing N would make the lines less thick but would not change their overall form, that is derived from the frequency characteristics of the stochasticity.

Conclusion

We have introduced and applied various advanced statistical model checking techniques to an osillatory biological system and have demonstrated their relevance to runtime verification. We have used two state of the art tools: Uppaalsmc and Plasma. Uppaal is a mature general-purpose platform based on timed automata, while Plasma is a relatively new tool that allows domain-specific languages to describe Markov chain models. Table 1 summarizes the aspects of both tools and shows that Plasma is targeted for domain specific languages and the final result analysis (such as frequency analysis) is done using external tools, whereas Uppaal provides generic features with integrated visualization of predetermined concepts. Although there are fundamental differences in the underlying semantics of the two tools, we have shown that they are united by the properties and problems of verification. In particular, frequency domain analysis and rare events are subjects of ongoing joint research.

Fig. 1 :Fig. 2 :

 12 Fig. 1: System of ordinary differential equations describing the genetic oscillator example.Fig. 2: Petri net representation of the initial state of the genetic oscillator example.

  where a, b ∈ IN, a < b and alpha is a proposition of the model. In the case of the genetic oscillator we consider in this paper, the propositions are numeric comparisons of the variables in question. An expression as φ 1 U[a; b]φ 2 means that φ 1 should be true until φ 2 is true and this should occur between a and b time units. The expression φ 1 R[a; b]φ 2 means φ 2 should be true until either b time units has passed or both φ 1 and φ 2 is true and that occurs between a and b time units. Xφ means that φ should be true in the next state of the system. The remaining operators have their standard interpretation from propositional calculus. The derived eventuality operator ♦ [a;b] φ is introduced as short for ⊤U [a;b] φ and the always operator [a;b] φ as short for ⊥R [a;b] φ.

Fig. 4 :

 4 Fig. 4: The bouncing ball and the hitting piston (a), the automata for (b) the ball and (c) the piston, and (d) 5 trajectories of the ball in function of time.

  gammaR=0, betaA=0, deltaA=0, ... , A=0, C=0, R=0, DR=1, ... gammaR'==0 && betaA'==0 && deltaA'==0 && ... && A'== betaA*MA+thetaA*D_A+thetaR*D_R -A*(gammaA*DA+gammaR*DR+gammaC*R+deltaA) && R'== betaR*MR-gammaC*A*R+deltaA*C-deltaR*R && C'== gammaC*A*R-deltaA*C && ...(a) ODE model.

Fig. 5 :

 5 Fig. 5: Snippets of the Uppaal models of the genetic oscilator.

k→

  C + D have the concrete form A + B k-> C + D with guarded command semantics (A > 0 ∧ B > 0, kAB, A = A -1; B = B -1; C = C + 1; D = D + 1).

  . With initial conditions A = 1000, B = 1000, C = D = E = 0, we then consider the property Pr[♦ D ≥ 470] that has a probability of approximately 2 × 10 -10

  Stochastic model simulation plot.

Fig. 7 :

 7 Fig. 7: Uppaal-smc simulations: simulate 1 [<=75] { A, C, R }. than it normally would, thus we ignore the first period in our measurements. The amplitude can be measured from the plot directly for the ODE model: A -1375, C -2183 and R -1717. The amplitude is not fixed for stochastic model, and thus the distribution of probable amplitude values is estimated instead. The start of the monitoring is constrained with a variable v gaining value 1 only after 15 time units, effectively ignoring the first peaks. The results of 2000 simulations of 75 time units are shown in Fig. 8. The average amplitudes are a bit larger than in ODE model.We can also estimate periods of various signals by measuring the time distance between peaks. For this purpose we add signal peak monitor-process gener-

Fig. 8 :

 8 Fig. 8: Estimated amplitude: E[<=75; 2000](max: A*v,C*v,R*v).

Fig. 9 :Fig. 10 :

 910 Fig. 9: Estimated periods: Pr[x<=100](<> secondPeak.ACCEPT).

Fig. 11 :

 11 Fig. 11: Average frequency spectra of protein R. δt = 2h, K = 12500. N = 100 (N = 1) for the stochastic (deterministic) model.

Fig. 12 :

 12 Fig. 12: Average frequency spectra of protein C. δt = 2h, K = 12500. N = 100 (N = 1) for the stochastic (deterministic) model.

Table 1 :

 1 Tool comparison summary.

	Aspect	Plasma	Uppaal
	Models	Domain specific	Stochastic hybrid automata
	Semantics DTMC/CTMC	Stochastic hybrid transition systems
	Properties MITL, rare events (Weighted) MITL
	Results	Exported data	Generic visualizations

i.e., a property whose semantics is trace-based.

⋆ Work partially supported by VKR Centre of Excellence -MT-LAB and by the CREATIVE project ESTASE.


Fabrizio Biondi

Axel Legay

Quantitative Anonymity Evaluation of Voting Protocols

In an election, it is imperative that the vote of the single voters remain anonymous and undisclosed. Alas, modern anonymity approaches acknowledge that there is an unavoidable leak of anonymity just by publishing data related to the secret, like the election's result. Information theory is applied to quantify this leak and ascertain that it remains below an acceptable threshold. We apply modern quantitative anonymity analysis techniques via the state-of-theart QUAIL tool to the voting scenario. We consider different voting typologies and establish which are more effective in protecting the voter's privacy. We further demonstrate the effectiveness of the protocols in protecting the privacy of the single voters, deriving an important desirable property of protocols depending on composite secrets.

Introduction

Voting is the backbone of the democratic process [START_REF] Norris | Electoral Engineering: Voting Rules and Political Behavior[END_REF]. To be effective, a voting system must allow the voters to freely express their opinion and elect the public officials that will represent them in the government. An effective voting system guarantees that each vote is counted exactly once, that no malicious agent can tamper with the results of the vote, and that no vote can be traced back to the voter who cast it.

Various traditional and electronic voting systems have been proposed to assure such guarantees. The use of cryptography and certification authorities can guarantee that only eligible voters can vote and that their vote is counted exactly once, and the production of fake credentials can safeguard voters against being coerced to reveal their vote [START_REF] Clarkson | Civitas: Toward a secure voting system[END_REF]. The anonymity of the vote is harder to guarantee; current proposals include assumptions on the absolute anonymity of the voting channels [START_REF] Clarkson | Civitas: Toward a secure voting system[END_REF][START_REF] Gritzali | Principles and requirements for a secure e-voting system[END_REF] or expect enterprises and universities to provide public proxy servers to hide the IP address of the voter [START_REF] Chen | The design of a secure anonymous internet voting system[END_REF]. The problem with these approaches to anonymity is that the anonymity of a vote is considered a qualitative, yes/no property, verifying whether it is possible for an attacker to infer any amount of information about the identity of the voters; this is known as the possibilistic approach. The probabilistic approach instead considers anonymity as a quantity that can be decreased by the attack of an external agent and by other factor in the voting process including the magnitude of the electoral seat, the electoral formula used, the results of the elections and the number of candidates. Since none of these factors completely compromise anonymity, a qualitative technique has to either ignore them or consider any voting protocol unsafe.

Current approaches to anonymity consider a secret, like the identity of the caster of a vote, as a quantitative amount of information, and use information theory to quantify how much of this secret information is inferred by a malicious attacker [START_REF] Chatzikokolakis | Anonymity protocols as noisy channels[END_REF][START_REF] Chen | Quantifying maximal loss of anonymity in protocols[END_REF]. This amount is called information leakage. The qualitative approach tags as insecure even a negligible amount of loss of information, in practice considering any real system insecure except under very strong assumptions. On the other hand, the quantitative approach allows the analyzer to determine a bound above which a loss of anonymity is considered noteworthy. Quantitative anonymity analysis has been applied among others to study the trade-off between anonymity and utility of operations on databases [START_REF] Alvim | Differential privacy: On the trade-off between utility and information leakage[END_REF] and to define information-theoretical bounds to differential privacy [START_REF] Barthe | Information-theoretic bounds for differentially private mechanisms[END_REF]. To the best of our knowledge, these techniques have not been applied to voting protocols.

Completely automated tools have been created to quantify information leakage for any secret-dependent protocol. Previously we have introduced a theoretical framework to model protocols with Markov chains and efficiently and precisely compute their leakage [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF]. We implemented the approach in the QUAIL tool, the first tool able to perform an arbitrary-precision leakage analysis of a non-deterministic secret-dependent protocol [START_REF] Biondi | QUAIL: A quantitative security analyzer for imperative code[END_REF].

In this work we will use QUAIL to analyze different typologies of electoral formulae. QUAIL considers the combined votes of the voters as a precise amount of secret bits and quantifies precisely how many of these bits are inferred by an attacker able to read the published results of the elections, i.e. the information leakage. Since these results are public, there is no way to avoid this loss of anonymity. Consequently any qualitative method that claims to perfectly guarantee anonymity of the voters is ignoring this nonnegative information leakage. By quantifying exactly this amount we establish a lower bound on the amount of anonymity that any implementation or formula can guarantee.

We study two very general typologies of electoral formulae: Single Preference formulae, where each voter expresses a preference for a candidate, and Preference Ranking formulae, where each voter ranks all candidates from the best to the worst. Our results are valid for any electoral formula in the typologies. This classification is traversal to the common division in proportional and majoritarian systems, as it depends only on the way the vote is expressed.

We consider that the secret is not a single entity, but a composite of the secret votes. If we have 10 secrets and each is 2 bits and the leakage of the system is 2 bits, we need to identify whether the leaked information corresponds to the secret of one of the voters or only to general information about the result of the vote. The theoretical bases for the study of leakage of composite secrets are very recent [START_REF] Alvim | Measuring information leakage using generalized gain functions[END_REF], and to the best of our knowledge no large scale study of practical cases has been published.

Our analysis shows the exact impact of the magnitude of the electoral seat and the number of candidates on the anonymity of the vote. The results suggest the data from seats of low magnitude, like hospital seats, should be aggregated before publication to protect the voters' anonymity. Also we show that the leakage on a single voter's secret is strictly less than the leakage on all secrets divided by the number of voters, proving that the voting protocols analyzed are effective in protecting the single voters' secrets. This is a generally desirable property for protocols on composite secrets.

The rest of the paper is organized as follows: Section 2 introduces common concepts in probability and information theory and Section 3 details the leakage theory and how it is implemented in the QUAIL analyzer. Section 4 presents the two typologies of electoral formulae we analyze and in Section 5 we discuss the results obtained. Section 6 explains which problems we are facing in obtaining more results and what steps we are taking to solve them.

Background

We define common concepts in probability theory and information theory that are used through the paper. We refer to books on the subject [START_REF] Cover | Elements of information theory[END_REF] for the basic definitions on probability theory. We call X a discrete random variable and X a discrete stochastic process, i.e. an indexed infinite sequence of discrete random variables (X 0 .X 1 , X 2 , ...) ranging over the same sample space S. The index of the random variables in a stochastic process can be understood as modeling a concept of discrete time, so X k is the random variable representing the system at time unit k.

Markov Chains

A discrete stochastic process is a Markov chain C = (C 0 , C 1 , C 2 , ...) iff ∀k ∈ N. P (C k |C k-1 , ..., C 0 ) = P (C k |C k-1
). A Markov chain on a sample space S can also be defined as follows: Definition 1. A tuple C = (S, s 0 , P ) is a Markov Chain (MC), if S is a finite set of states, s 0 ∈ S is the initial state and P is an |S| × |S| probability transition matrix, so ∀s, t ∈ S. P s,t ≥ 0 and ∀s ∈ S. t∈S P s,t = 1.

We call π (k) the probability distribution vector over S at time k and π (k) s the probability π k (s) of visiting the state s at time k. This means that, considering a Markov chain C as a time-indexed discrete stochastic process (C 0 , C 1 , ...), we write π (k) for the probability distribution over the random variable C k . Since we assume that the chain starts in state s 0 , then π (0) s is 1 if s = s 0 and 0 otherwise. Note that π (k) = π 0 P k , where P k is matrix P elevated to power k, and P 0 is the identity matrix of size |S| × |S|.

A state s ∈ S is absorbing if P s,s = 1. In the figures we will not draw the looping transition of the absorbing states, to reduce clutter. We say that a Markov chain is one-step if all states except the starting state s 0 are absorbing. We will usually refer to one-step Markov chains as C.

Let ξ(s, t) denote the expected residence time in a state t in an execution starting from state s given by ξ(s, t) = ∞ n=0 P n s,t . We will write ξ s for ξ(s 0 , s). We will enrich our Markov chains with a finite set V of natural-valued variables, and for simplicity we assume that there is a very large finite bit-size M such that a variable is at most M bit long. We define an assignment function A : S → [0, 2 M -1] |V| assigning to each state the values of the variables in that state. We will use the expression v(s) to denote the value of the variable v ∈ V in the state s ∈ S.

Given a Markov chain C = (S, s 0 , P ) let a discrimination relation R be an equivalence relation over S. We use discrimination relation to quotient one-step Markov chains: Definition 2. Given a one-step Markov chain C = (S, s 0 , P ) and a discrimination relation R over S, we define the quotient C/R of C over R as the one-step Markov chain C/R = ( S, s0 , P ) where -S is the set of equivalence classes of S induced by R; -s0 is the equivalence class of s 0 ; for each equivalence class s ∈ S , Ps0,s = s∈s P s0,s and Ps,s = 1.

Information Theory

The entropy of a probability distribution is a measure of the unpredictability of the events considered in the distribution [START_REF] Shannon | A mathematical theory of communication[END_REF]. Definition 3. [START_REF] Cover | Elements of information theory[END_REF] Let X and Y be two random variables with probability mass functions p(x) and p(y) respectively and joint pmf p(x, y). Then we define the following non-negative real-valued functions:

-Entropy H(X) = -x∈X p(x) log 2 p(x) -Joint entropy H(X, Y ) = -x∈X y∈Y p(x, y) log 2 p(x, y) -Conditional entropy H(X|Y ) = -x∈X y∈Y p(x, y) log 2 p(x|y) = = y∈Y p(y)H(X|Y = y) = y∈Y p(y) x∈X p(x|y) log 2 p(x|y) = = H(X, Y ) -H(Y ) (chain rule) -Mutual information I(X; Y ) = x∈X y∈Y p(x, y) log 2 p(x,y) p(x)p(y) = = H(X) + H(Y ) -H(X, Y ) ≤ min(H(X), H(Y ))
Since every state s in a MC C has a discrete probability distribution over the successor states we can calculate the entropy of this distribution. We will call it local entropy, L(s), of s: L(s) = -t∈S P s,t log 2 P s,t . Note that L(s) ≤ log 2 (|S|).

As a MC C can be seen as a discrete probability distribution over all of its possible traces, we can assign a single entropy value H(C) to it. The global entropy H(C) of C can be computed by considering the local entropy L(s) as the expected reward of a state s and then computing the expected total reward of the chain [START_REF] Biondi | Maximizing entropy over Markov processes[END_REF]: H(C) = s∈S L(s)ξ s If a Markov chain is one-step its entropy corresponds to the local entropy of the initial state s 0 .

Information Leakage of Markov Chains

Theoretical background

We use information theory to compute the amount of bits of a secret variable h that can be inferred by an attacker able to observe the value of an observable variable o after the termination of a protocol. We call this amount Shannon leakage or just leakage, and it corresponds to the mutual information between the distribution on the secret and the distribution on the observable variable. This analysis assumes the worst possible attacker: the attacker has access to the source code of the protocol and to unlimited computational power.

We will model system-attacker scenarios with Markov chains in which to each state we associate a unique assignment of values to all variables [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF]. Then we define leakage as follows: Note that to compute leakage we need to have a prior probability distribution over the secret, modeling what the attacker knows before observing the observable output of the protocol. We will assume for simplicity that the attacker knows the possible values of the secret, since he can read the source code and verify which kind of variable holds it, but has no additional information about it.

The modeling of a system-attacker scenario as a Markov chain starts by dividing the system's variables in private and public variables. Private variables, including the secret variable h, are the ones whose value is not defined at compilation time. In each state of the Markov chain a set of allowed values is assigned to each private variable. Public variables, including the observable variable o and the program counter pc, are variable whose value is known to the analyst. On each state a given value is assigned to each public variable.

Given the source code of the system and a prior distribution over the private variables, we have enough information to build the Markov chain semantics C of the protocol, since for each state we can determine its successor states and the corresponding transition probabilities. We show a simple example, and refer to [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF] for the complete semantics. Let secret variable h be a secret bit, observable variable o an observable bit and public variable r a random bit being assigned the value 0 with probability 0.75 and 1 otherwise. We assign to o the result of the exclusive OR between h and r and terminate. We want to quantify the amount of information about h that can be inferred by knowing the value of o. The Markov chain semantics C for the example is shown in Fig. 1a. Each state is enriched with information about the allowed values of private variables and the values of public variables, e.g. in state S 1 secret variable h can be either 0 or 1 and public variable r is 0.

Subsequently we need to model the fact that the attacker has to wait for the protocol to terminate to read the observable output. Equivalently, we can say that the attacker is not aware of the internal behavior of the system; this is modeled by hiding in the Markov chain model the internal states, i.e. all states except the initial state and the absorbing states. We call the resulting Markov chain the observable reduction C. The observable reduction for the example is shown in Fig. 1b.

Note that the observable reduction is a one-step Markov chain, so we can compute quotients on it following Definition 2. To compute the leakage we need to compute three quotients from the observable reduction: Finally we compute leakage as the mutual information I(O; h), as between the secret and observable variable, as explained in Definition 4. To compute it we apply the formula

I(X; Y ) = H(X) + H(Y ) -H(X, Y ) from Definition 3, obtaining I(o; h) = I(C |o ; C |h ) = H(C |o ) + H(C |h ) -H(C |(o,h) )
meaning that the leakage can be computed as the sum of the entropies of the secret and observable quotient minus the entropy of the joint quotient. In our example we have 1 andH(C |(o,h) ) = 1.8112..., so we conclude that the program leaks 1 + 1 -1.8112... ≈ 0.1887 bits, or 18.87% of the secret.

(C |o ) = 1, H(C |h ) =

QUAIL Implementation

The QUAIL tool quantifies the Shannon leakage of a probabilistic process in a fully automated way. The user just provides the source code for the process in the QUAIL imperative language, specifying the size of the variables and whether each variable is public, private, secret or observable. QUAIL computes the expected amount of information about the secret variables that an attacker is able to infer by knowing the values of the observable variables after the protocol has terminated and produced output, implementing the theory presented in Section 3.1.

The QUAIL language is a WHILE language enriched with for loops, multidimensional arrays and constant declarations. For simplicity all variables are integer variables; the bit length of each variable is defined by the coder at declaration time. Also, assignment to private or secret variables is not allowed, and all variables have to be declared at the beginning of the program. We refer to [1] for the source code and full semantics of QUAIL.

QUAIL is able to handle secret and private variables with a large number of possible values. Alas, the representation of arbitrary probability distributions on such variables on a real machine is untreatable: a probability distribution over a 64-bit variable is composed of 2 64 ≈ 1/8 • 10 19 rational numbers. For this reason QUAIL does not allow the user to define arbitrary prior distributions over the secret, and always assumes that the prior distribution is uniform. The QUAIL source code for the bit XOR example presented in Section 3.1 is shown in Fig. 2. In lines 1 to 3 the variables are declared, then in line 4 variable r is assigned value 0 with probability 0.75 and 1 otherwise, and in line 5 variable o is assigned with the exclusive OR of variables h and r.

We save the source code in Fig. 2 in a file bitxor.quail and invoke QUAIL with the command ./quail bitxor.quail -v 0 -p 5 where -v 0 suppresses all output except for the leakage result and -p 5 specifies that we want 5 significant digits in the answer. QUAIL outputs 0.18872 showing that the program leaks ≈ 0.18872 bits of the secret, in accordance to what we computed theoretically in Section 3.1.

Modeling Voting Protocols

In an election, each voter is called to express his preference for the competing candidates. The voting system defines the way the voters express their preference: either on paper in a traditional election, or electronically in e-voting. The voting system also comprehends the additional procedures enforced to guarantee that the voters can vote freely, that they can verify that their vote has been counted and that their vote remains confidential.

After the votes have been cast, the results of the vote are published, usually in an aggregated form to protect the anonymity of the voters. Finally, the winning candidate or candidates is chosen according to a given electoral formula.

In this section we present two different typologies of voting, representing two different ways in which the voters can express their preference: in the Single Preference protocol the voters declare their preference for exactly one of the candidates, while in the Preference Ranking protocol each voter ranks the candidate from his most favorite to his least favorite. Since each protocol we model is concerned only about how the votes are expressed and counted and what results are published, each protocol models a number of electoral formulae. For the same reason, the models are valid both for uninominal and multinominal elections.

Single Preference

The Single Preference protocol typology models all electoral formulae in which each of the N voters expresses one vote for one of the C candidates, including plurality and majority voting systems and single non-transferable vote [START_REF] Norris | Electoral Engineering: Voting Rules and Political Behavior[END_REF]. The votes for each candidate are summed up and only the results are published, thus hiding information about which voter voted for which candidate. The candidate or candidates to be elected are decided according to the electoral formula used.

Secret and observable encoding. The secrets and observables are modeled by the following lines of QUAIL code: The secret is an array of integers with a value for each of the N voters. Each value is a number from 0 to C-1, representing a vote for one of the C candidates. The observable is an array of integers with the votes obtained by each of the C candidates. The full model for this protocol is shown in the Appendix due to space constraints.

Formal leakage computation. The protocol is simple, and its information leakage can be computed formally, as shown by the following lemma: Lemma 1. The information leakage for the Single Preference protocol with n voters and c candidates corresponds to

log 2 c n - 1 c n k1+k2+...+kc=n n k 1 + k 2 + ... + k c log 2 n k 1 + k 2 + ... + k c
The proof for Lemma 1 is in the Appendix. While the lemma characterizes the solution computed by QUAIL for this case, it is very hard to find such a characterization for any process, so in general QUAIL is the best way to obtain a result. We run QUAIL with the command ./quail single_preference.quail -v 0 -p 5 with the same parameters we used in Section 3.2 and obtain 1.8112 showing that the leakage of the Single Preference protocol for 3 voters and 2 candidates is ≈ 1.8112 bits.

Preference Ranking

The Preference Ranking protocol typology models all electoral formulae in which each of the n voters expresses an order of preference of the c candidates, including the alternative vote and single transferable vote systems [START_REF] Norris | Electoral Engineering: Voting Rules and Political Behavior[END_REF]. In the Preferential Voting protocol the voter does not express a single vote, but rather a ranking of the candidates; thus if the candidates are A, B, C and D the voter could express the fact that he prefers B, then D, then C and finally A. Then each candidate gets c points for each time he appears as first choice, c -1 points for each time he appears as second choice, and so on. The points of each candidate are summed up and the results are published.

Secret and observable encoding. The secrets and observables are modeled by the following lines of QUAIL code: The secret is an array of integers with a value for each of the N voters. Each value is a number from 0 to C!-1, representing one of the possible C! rankings of the C candidates. The observable is an array of integers with the points obtained by each of the C candidates. The full model for this protocol is shown in the Appendix due to space constraints.

Experimental Results

We discuss some of the initial results we have obtained by analyzing the Single Preference and Preference Ranking voting protocols. Discussion. Since the posterior entropy measures how hard it would be for an attacker to learn the secret after observing the results of the voting, we focus on it as the measure of how confidential the votes are after the attack. Note that posterior entropy increases less than linearly with the number of targets, so if we want to learn both the votes of voters Alice and Bob it is more convenient to consider the two votes as a single composite secret than to try to learn the two votes separately. Note also that the posterior entropy when all voters are targets is the same as the posterior entropy when all voters except one are. This is because if all votes except one are known, the last one can be inferred immediately by checking the election's results. It should also be noted that the percentage of information leaked I O,h % increases with the number of targets, again sublinearly. We argue that this is a desirable property in a protocol designed to protect a secret composed of several subsecrets with the same importance. The property ensures that it is not more convenient for the attacker to try to infer separately every single secret instead of the whole composed secret, so even if it takes less time to infer the secret of 1 target out of 3, the time it would take to infer all 3 secrets one by one is larger than the time required to infer them all at the same time, as we explained in the paragraph above. This guarantees that the posterior entropies for multiple targets are in fact sound. This property also forces the attacker to decide beforehand exactly how many votes he needs to discover to minimize the time needed for the attack.

3V-2C HO H

h H O,h I O,h I O,h % pH h 1 1.

Challenges

Problem size

The examples we analyzed consider only a small number of voters and candidates. The algorithm for the precise computation of information leakage is exponential in the size of the secret and the size of the secret grows with the number of voters and candidates, thus QUAIL and any other tool are too slow to analyze large cases. Analyzing the Preference Ranking protocol also requires more time than analyzing the Single Preference protocol, since the former protocol is more complex than the latter and has a larger secret size.

To solve this problem we are implementing a statistical analyzer in QUAIL able to simulate the execution of the protocols a large number of times and to approximate the information leakage value by analyzing the collected data. The statistical analyzer is able to solve larger problems than the standard QUAIL algorithm, since it does not have to analyze the whole space of possible program executions.

g-leakage

Information leakage quantifies the loss of information on the whole secret. We have analyzed the behavior of leakage when we consider only some of the votes to be the secret we are interested in, showing how the protocol is efficient in hiding the secret of the single voters.

Recently an extension of information leakage, called g-leakage [START_REF] Alvim | Measuring information leakage using generalized gain functions[END_REF], has been proposed exactly to deal with cases in which different subsets of the secret bits have different values for the attacker, as is the case with composite secrets. g-leakage is very general, since it allows for any gain function to be used in top of leakage computation; for this reason it has not yet been implemented in any leakage analysis tool.

We are working to extend QUAIL with g-leakage computation capabilities, allowing us to encode more naturally problems with composite secrets like voting protocols. We expect that the results on the efficiency of the protocols in protecting the single votes will be coherent with the results presented in this paper.

Implementation Analysis

The protocols we analyze model the abstract behavior of voting systems, giving us theoretical lower bounds on the amount of information leaked by publishing the results of the elections. It would be interesting to compare these theoretical results with actual implementations of voting systems, to evaluate how effective the real systems are in guaranteeing anonymity. Off-the-shelf systems are obviously not written in QUAIL language, so a tool capable of analyzing C or Java code like the ones developed by Phan and Malacaria [START_REF] Phan | Abstract model counting: a novel approach for quantification of information leaks[END_REF] or by Chothia et al. [START_REF] Chothia | A tool for estimating information leakage[END_REF] would have to be used.

A Appendix

Proof (of Lemma 1). Call h the composite secret about the votes of all voters and O the observable output of the system. Remember that information leakage corresponds to the difference between the prior entropy on the secret H(h) and the posterior entropy on the secret after observing the observable output of the system H(h|O) = Preference Ranking The model for the Preference Ranking model is shown on Fig. 3 on the right. Constant N represents the number of voters and constant C the number of candidates. The observable variable result is an array with the total points obtained by each candidate, and the secret variable vote is an array with the preference ranking of each voter. For each voter the secret has c! possible different votes, corresponding to the possible complete orderings of the c candidates. The secret vote of each voter is encoded as a number from 0 to c! -1, and then transformed in a preferential order with the algorithm in lines 30-57. The points for each candidate are counted by summing the points given by the preference ranking of each voter.

Fig. 1 :

 1 Fig. 1: Bit XOR example: a) Markov chain semantics C. b) Observable reduction C. c) Joint quotient C |(o,h) . d) Secret's quotient C |h . e) Observer's quotient C |o .

  Joint quotient The joint quotient process C |(o,h) models the joint behavior of the secret and observable variables. It is shown in Fig. 1c. Secret's quotient The secret's quotient process C |h models the behavior of the secret variable. It is shown in Fig. 1d. Observer's quotient The observer's quotient process C |o models the behavior of the observable variable. It is shown in Fig. 1e.

1 secret int1 h; 2 Fig. 2 :

 22 Fig. 2: Bit XOR example: source code.

  secret array [ N ] of i n t 3 2 v o t e : = [ 0 , C-1]; observable array [ C ] of i n t 3 2 r e s u l t ;

  secret array [ N ] of i n t 3 2 v o t e : = [ 0 ,C! -1]; observable array [ C ] of i n t 3 2 r e s u l t ;

  o∈O P (o)H(h|O = o).The secret h has c n possible values, thus can be encoded in log 2 c n bits, thus the prior entropy H(h) corresponds to log 2 c n .For the posterior entropy H(h|O), consider that the possible votes on the candidates form a multinomial distribution, thus the probability P (o) of a given outcome o ∈ O is 1 /c n and the conditional posterior entropy H(h|O = o) where k i is the amount of votes to candidate 1 ≤ i ≤ c is n k1+k2+...+kc . We conclude thatI(O, h) = H(h) -o∈O P (o)H(h|O = o) = log 2 c n -1 c n k1+k2+...+kc=n n k 1 + k 2 + ... + k c log 2 n k 1 + k 2 + ... + k cVoting Protocol Models Single Preference The model for the Single Preference protocol is shown on Fig.3on the left. Constant N represents the number of voters and constant C the number of candidates. The observable variable result is an array with the total votes expressed for each candidate, and the secret variable vote is the array with the preference of each voter. The rest of the code just sums up the votes for each candidate.

Table 3 :

 3 Single Preference voting protocol: entropies and leakage on varying the number of voters, candidates and target voters.Analysis of Single Preference with variable number of targets. In Table3we give detailed results of the analysis of the Single Preference voting protocol, on varying the number of voters, candidates and target voters. The code in bold on the top left corner of a table shows how many voters and candidates are being considered in the experiment, e.g. 4V-2C means 4 voters and 2 candidates. The column of the left represents the number of target voters for the experiment, i.e. how many votes is the attacker trying to infer. The table reports the following values:H O is the entropy of the observer's quotient H h is the (prior) entropy of the secret's quotient H O,h is the entropy of the joint quotientI O,h = H O + H h -H O,h is the information leakage I O,h % = I O,h /H his the percentage of the secret that has been leaked pH h = H h -I O,h is the (expected) posterior entropy of the secret, i.e. the amount of secret that has not been leaked

			3V-3C HO H h H O,h I O,h I O,h % pH h
		81 1.00 2.50 0.31 31.1% 0.69	1	3.12 1.58 4.08 0.62 39.1% 0.96
	2	1.81 2.00 3.00 0.81 40.5% 1.19	2	3.12 3.16 4.75 1.53 48.5% 1.63
	3	1.81 3.00 3.00 1.81 60.3% 1.19	3	3.12 4.75 4.75 3.12 65.6% 1.63
	4V-2C HO H h H O,h I O,h I O,h % pH h	4V-3C HO H h H O,h I O,h I O,h % pH h
	1	2.03 1.00 2.81 0.21 21.0% 0.79	1	3.57 1.58 4.70 0.45 28.7% 1.13
	2	2.03 2.00 3.50 0.53 26.5% 1.47	2	3.57 3.16 5.67 1.07 34.0% 2.09
	3	2.03 3.00 4.00 1.03 34.0% 1.97	3	3.57 4.75 6.33 1.99 41.9% 2.76
	4	2.03 4.00 4.00 2.03 50.0% 1.97	4	3.57 6.33 6.33 3.57 56.3% 2.76
	5V-2C HO H h H O,h I O,h I O,h % pH h	5V-3C HO H h H O,h I O,h I O,h % pH h
	1	2.19 1.00 3.03 0.16 16.7% 0.84	1	3.93 1.58 5.16 0.35 22.3% 1.23
	2	2.19 2.00 3.81 0.38 19.3% 1.62	2	3.93 3.16 6.29 0.80 25.5% 2.36
	3	2.19 3.00 4.50 0.69 23.2% 2.31	3	3.93 4.75 7.25 1.43 30.1% 3.32
	4	2.19 4.00 5.00 1.19 30.0% 2.81	4	3.93 6.33 7.92 2.34 37.0% 3.99
	5	2.19 5.00 5.00 2.19 44.0% 2.81	5	3.93 7.92 7.92 3.93 49.6% 3.99