1 Introduction
Biological systems at the single cell level are inherently stochastic. Molecules inside the cells perform random movements (random walk) and the reactions among them may occur when the probability of collision is high enough. The number of molecules of each species at each time point is therefore a random process: assuming instantaneous reactions, this process can be modelled as a Markovian (i.e. memoriless) discrete state, continuous time process. When the number of molecules of each species involved is large, so that many reactions happen in any small interval of time, stochastic effects can be neglected. However, if the concentration of the molecules (of at least some of the species) is low the stochasticity plays an important role and must be taken into account. For this reason, stochastic models such as ContinuousTime Markov Chains (CTMC) [25] and Stochastic Hybrid Automata [15] are particularly powerful and suitable formalisms to model and to reason about biological systems defined as stochastic systems over time.
A classical question in formal modelling is to calculate the probability that a behaviour, expressed in terms of a certain temporal logic formula, may occur in a given stochastic process, with specified parameters. Probabilistic Model Checking [4, 5] (PMC) is a wellestablished verification technique that provides a quantitative answer to such a question. The algorithm used to calculate this probability [35] produces the exact solution, as it operates directly on the structure of the Markov chain. Despite the success and the importance of PMC, this technique suffers some computational limitations, either due to state space explosion or to the difficulty (impossibility) in checking analytically formulae in specific logics, like Metric Temporal Logic (MTL) [5, 17]. Furthermore, PMC provides only a quantitative measure of the satisfability (yes/no answer) of a temporal logic specification (i.e., the probability of the property being true).
However, especially when we deal with stochastic models, the notion of satisfability may be not enough to determine the capacity of a system to mantain a particular emergent behaviour unaffected by the uncertainty of the perturbations due to its stochastic nature or by possible small changes in the model parameters. A similar issue also arises when considering the satisfability of a property by deterministic dynamical systems which may be subject to extrinsic noise or uncertainty in the parameter. To address this question in the deterministic case, researchers from the verification community have proposed several notions of temporal logic based robustness [24, 28, 42], providing suitable definitions of distance between a trajectory of a system and the behavioural property of interest, expressed in terms of a temporal logic formula. These effectively endow the logic of interest with quantitative semantics, allowing us to capture not only whether a property is satisfied but also how much it is satisfied. A similar notion of robustness for stochastic models would clearly be desirable but, to our knowledge, has not been formalised yet.
The contributions of this paper are twofold. First,we provide a simulationbased method to define a notion of robust satisfability in stochastic models. Simulationbased approaches, such as statistical model checking [44]
, can be be used to estimate for a stochastic model the
robust satisfability distribution for a given temporal logic formula, with a guarantee of asymptotic correctness. This distribution is the key to understand how the behaviour specified by the logic temporal formula is unaffected by the stochasticity of the system. In particular, in this paper we consider two important indicators of this distribution: the robustness average and the conditional robustness average on a formula being true or false. We discuss how to compute the robust satisfability distribution and its indicators on two biological examples. Second, we show how to combine these indicators with the satisfaction probability to address the system design problem, where the goal is to optimize (few) control parameters of a stochastic model in order to best maximize these three indicators. The proposed approach takes advantage of Gaussian Process Upper Confidence Bound (GPUCB) algorithm introduced in [43].The paper is structured as follows: in Section 2 we introduce the background material. In Section 3 we discuss the robustness of stochastic models using the quantitative semantics of the Signal Temporal Logic (STL). In Section 4 we present some experimental results for the robustness of STL formulae for two stochastic models that we have chosen as our case studies: the Schlögl system and the Reprissilator. In Section 5 we show an application of the robust semantics to the system design problem. The related works and the final discussion are in Section 6.
2 Background
Markov Population Models.
The simplest class of stochastic processes we will consider are Continuous Time Markov Chains (CTMC) [25] that describe population processes (PCTMC). A population process intuitively is a system in which agents or objects of different kinds, and with different internal states, interact together. The classical example are biochemical and genetic networks, but other population processes include ecological systems, computer networks, and social systems.
We will describe PCTMC by the simple formalism of biochemical reaction networks. The state of the system is described by a vector
of integervalued variables , each counting the number of entities of a given class or species. The dynamics of this system is specified by a set of reactions , which can be seen as description of events changing the state of the system. Each reaction is of the formwhere is a reactant and is a product (they are both variables of ), and , are the stoichiometric coefficients, i.e. the amount of agents/ entities consumed or produced by the reaction. Stoichiometric information of a reaction can be condensed into an update vector , giving the net change in population variables due to : , where equals one in position and zero elsewhere. Additionally, each reaction has an associated rate function giving the rate of the transition as a function of the global state of the system.
From a set of reactions and species , we can easily derive the formal representation of a CTMC in terms of its infinitesimal generator matrix, see for instance [10]. Here we just recall that the state space of the CTMC is (or a proper subset, if any conservation law is in force). Such CTMC can be simulated with standard algorithms, like SSA [30].
Fluid Approximation
From a Markov population model, we can easily construct an alternative semantics in terms of Ordinary Differential Equations (ODE), assuming variables
to be continuous and interpreting each rate as a flow, thus obtaining the vector fielddefining the ODE . This equation, known as fluid approximation, can be shown to be a first order approximation of the average of the CTMC, and, under a suitable rescaling of the variables (dividing by the system size, which for biochemical reactions is just the volume), one can prove convergence of the CTMC to the solution of the ODE (see [10]) as populations and system size go to infinity. Intuitively, this ODE is a good description of the system behaviour when populations are large.
Stochastic Hybrid Automata
In many situations, it is not the case that all entities/ species in the model are present in large quantities. In such scenarios, fluid approximation can give poor results, yet dealing with CTMCs can be computationally unfeasible. An example are genetic regulatory networks, in which genes are modelled explicitly as a finite state machine [13]. In these cases, a better strategy is that of approximating continuously only some variables, keeping discrete the others. This reflects in the dynamics: some reactions will be converted into flows (generally those modifying only continuous variables), while the others will remain stochastic discrete events. This gives rise to a model that can be expressed in terms of a class of Stochastic Hybrid Automata (SHA, [13]) known as PiecewiseDeterministic Markov Processes [18]. Alternatives assuming a stochastic continuous dynamics have also been considered [38, 39]. More specifically, the SHA so obtained have discrete modes identified by the value of discrete variables. In between discrete transitions, the system evolves following the solution of the differential equation, whose vector field is modedependent (via the value of discrete variables). Discrete jumps happen at exponentially random distributed times, at a nonconstant rate that can depend on the continuous variables. After each jump, the value of discrete variables can change. Also continuous variables can be updated, even if we do not consider this possibility in this paper, see [13] for further details. Similarly to the fluid approximation case, we can see SHA models as the limit of CTMC, taking to the limit only the populations corresponding to continuous variables (under a suitable scaling of rates, see [9] for further details).
Signal Temporal Logic.
Temporal logic [40] provides a very elegant framework to specify in a compact and formal way an emergent behaviour in terms of timedependent events. Among the myriads of temporal logic extensions available, Signal Temporal Logic [37] (STL) is very suitable to characterize behavioural patterns in time series of real values generated during the simulation of a dynamical system. STL extends the densetime semantics of Metric Interval Temporal Logic [2] (MITL), with a set of parametrized numerical predicates playing the role of atomic propositions. STL provides two different semantics: a boolean semantics that returns yes/no depending if the observed trace satisfies or not the STL specification and a quantitative semantics that also returns measure of robustness of the specification. Recently, Donzé et. al [23] proposed a very efficient monitoring algorithm for STL robustness, now implemented in the Breach [20] tool. The combination of robustness and sensitivitybased analysis of STL formulae have been successfully applied in several domains ranging from analog circuits [32] to systems biology [21, 22], to study the parameter space and also to refine the uncertainty of the parameter sets. In the following we recall [24] the syntax and the quantitative semantics of STL that will be used in the rest of the paper. The boolean semantics can be inferred using the sign of the quantitative result (positive for true and negative for false).
Definition 1 (STL syntax)
The syntax of the STL is given by
where is a true formula, conjunction and negation are the standard boolean connectives, is a densetime interval with and is the until operator.
The atomic predicate is defined as , where , , , is the primary signal, and is a realvalued function known as the secondary signal.
The (bounded) until operator requires to hold from now until, in a time between and time units, becomes true. The eventually operator and the always operator can be defined as usual: ,
Definition 2 (STL Quantitative Semantics for space robustness)
where is the quantitative satisfaction function, returning a real number quantifying the degree of satisfaction of the property by the signal at time . Moreover, .
We stress here that the choice of the secondary signals is an integral part of the definition of the STL formula expressing the behaviour of interest. Different choices of secondary signals result in different formulae, hence in different robustness measures. We also remark that the robustness score of Definition 2 has to be interpreted as a weight of how much a given model (with fixed initial conditions and parameters) satisfies an STL behaviour. More precisely, its absolute value can be seen as a distance of the signal under consideration from the set of trajectories satisfying/ dissatisfying the formula [28], in STL trajectory are projected with respect to secondary signals [24]. In this sense, this measure is different from the more common sensitivitybased notions of robustness, like those discussed in [34], measuring the size of a region in the parameter space in which the system behaviour is roughly constant. However, sensitivity analysis and its related techniques can be applied to the robustness score of Definition 2.
3 Robustness of Stochastic Models
Consider a STL formula , with predicates interpreted over state variables of a PCTMC model . The boolean semantics of is readily extended to stochastic models as customary, by measuring the probability of the set of trajectories of the CTMC that satisfy the formula:
The rationale behind such definition is that a PCTMC model defines a probability distribution on the space of trajectories, which is usually obtained by applying the cylindric construction
[5]. Furthermore, the set of trajectories that satisfy/ falsify a formula is a measurable set, so that we can safely talk about its probability. In the following, we will refer to the space of trajectories as , and interpret the PCTMC modelas a random variable
over . In order to extend this definition to the robustness score, it is convenient to think of the set of trajectories that satisfy as a measurable function , such that if and only if . Then, we can define the random variable on induced by the PCTMC via as the Bernoulli random variable which is equal to 1 with probability . We can equivalently write:We can extend the robustness score to PCTMC models in the same way: given a trajectory , we can compute its robustness score according to Def. 2 and interpret as a function from the trajectories in to . This function is easily seen to be measurable (with respect to the algebra induced from the Skorokhod topology in ), and so it induces a realvalued random variable with probability distribution given by
Staten otherwise, if we apply the definition of robustness to a stochastic model, we obtain a distribution of robustness degrees. This distribution tells us much more than the standard probabilistic semantics, because it tells us “how much” a formula is true.
In particular, in this paper we will be interested in some statistics of this distribution, specifically the average robustness degree, and the average robustness conditional on a formula being true or false. The first quantity gives a measure of how strongly a formula is satisfied on average. The larger this number, the more robust is satisfaction. Most of the times, this number will be correlated with the satisfaction probability, yet we can have a large average satisfaction score even for a small probability of satisfaction. Better indicators of the intensity of satisfaction and dissatisfaction are the conditional averages, and . These are related to the average by the equation
which holds provided is zero.
One goal of this paper is to investigate to what extent these three synthetic indices are good descriptors of the robustness distribution, and how they can be exploited to do parameter synthesis for PCTMC models.
Reaction  rate constant  init pop 

4 Case Studies
In this section, we investigate experimentally the notion of robust semantics of STL formulae for stochastic models. We will consider two systems: the Schlögl system [31], a simple network of biochemical reactions exhibiting a bistable behaviour, and the Repressilator [26], a synthetic biological clock implemented as a network of gene regulations. More specifically, we consider a CTMC model of the Schlög system and a hybrid model of the Repressilator [11, 12], in order to illustrate the general applicability of the stochastic robust semantics introduced in Section 3.
4.1 Schlögl system
The Schlögl model is a simple biochemical network with four reactions, listed in Table 1. The rates of the reactions are computed according to the mass action principle for stochastic models [30]. Species and are considered to be present in large quantities, hence assumed constant. The characteristic of this system is to have, for certain parameter values, like the one shown in Table 1, a bistable behaviour. More specifically, the reaction rate ODE system has two stable steady states, and for this model the trajectories of the stochastic system starting from a fixed initial state can end up in one attractor or the other. The probability of choosing one stable state or the other depends on the position of relatively to the basin of attraction of the two equilibria. If we start close to its boundary, the bistable behaviour becomes evident, see Figure 1(a).
We now consider the property of eventually ending up in one basin of attraction, and express it with the STL formula
(4.1) 
stating that the system, after at most units of time, stabilises to a value which remains above for as long as units of time. In this formula, the predicate corresponds to the linear secondary signal . As can be seen from Figure 1(a), if the model is in the large equilibrium, then this property will be true, and false in the other case.
If we estimate the probability of the formula statistically, then we obtain the value (10000 runs, error at confidence level). However, this raw number does not tell us anything specific about bistability. A system stabilising just above the threshold 300, such that roughly of its trajectories cross it “frequently”, may satisfy the same formula with the same probability. However, the bimodal behaviour becomes evident if we look at the distribution of the robustness degree of the formula, see Figure 1(b) . Hence, the robustness score carries an additional amount of information with respect to the satisfaction probability of a STL formula. We stress that we are not comparing the robustness degree with the probability distribution of the CTMC : both the satisfaction probability of and its robustness are (unidimensional) quantities derived from , which are easier to compute and visualise.
In Figure 2, we investigate the behaviour of the average robustness degree, and its relationship with the satisfaction probability. In order to do this, we varied the threshold level in the formula (Figure 2(a)), and the rate constant (Figure 2(b)), and estimated statistically the average robustness degree and the satisfaction probability from 10000 runs for each parameter combination. As we can see these two quantities are correlated. When we vary the threshold, the correlation between satisfaction probability and robustness score is around 0.8386, while the dependency seems to follow a sigmoid shaped curve. In the second case, instead, the correlation between satisfaction probability and average robustness degree is 0.9718, with an evident linear trend.
Finally, we consider the conditional robustness degrees. For model parameters as in Table 1, the average robustness conditional on Formula (4.1) being true is 169.89, while the robustness conditional on the formula being false is 239.52 (see also Figure 1(b)). These two indicators estimate how robustly the system remains in the basin of attraction of each steady state.
4.2 Repressilator
The second case study is a genuine stochastic hybrid model of the Repressilator [26], a synthetic genetic clock composed of three genes expressing three transcription factors repressing each other in a cyclical fashion (see Figure 3(a)). The stochastic hybrid model we consider is taken from [11, 12]. In the model, we lump the transcription and translation in a single event, and model production and degradation of the protein as continuous flows. The binding and unbinding of transcription factors from gene promoters, instead, are modelled as discrete and stochastic events. As we can see in Figure 3(b), the model exhibit sustained oscillations, albeit with an irregular period. This happens for parameters giving a strong repression via a low unbinding rate.
In order to check for the presence of oscillations, we use the STL formula
(4.2) 
expressing the fact that low values of alternate to high values, with a period between and . The secondary signals are , , and so on. Here can be one of the three proteins of the Repressilator. In the next discussion, we focus on .
Again, the robustness score gives us a measure of the satisfaction/dissatisfaction of the formula. As we can see from Figure 4(a), the robustness degree shows a bimodal behaviour also in this case. In particular, in case the formula is false, it gives some degree of information on the amplitude of oscillations, and on the stability of the period (relatively to the formula parameters). In fact, a robustness value of, say, 50 can be obtained for instance if from a point in which , the system remains below for a whole (half) period of the oscillation (which is constrained to be in ). This can happen due to low amplitude or irregular period. In Figure 4(b), we plot the average robustness degree against the satisfaction probability, varying the property parameter , showing once again the correlation between the two quantities.
5 System Design
We now discuss an application of the robust semantics to the system design problem. The problem we want to tackle is the following:
given a population (hybrid) model, depending on a set of parameters , and a specification given by a STL formula, find the parameter combination such that system satisfied with probability at least as robustly as possible.
We will tackle this problem by:

rephrasing it as an unconstrained optimisation problem, where we seek to optimise the average robustness, using penalty terms to encode for probability constraints. More specifically, assuming we want to enforce the satisfaction probability to be at least , we add a penalty term of the form , if , and 0 otherwise, where controls the penalty intensity.

evaluating the function to optimise using statistical model checking with a fixed number of runs, usually set to 100;

solving the optimisation problem using an optimisation strategy for reinforcement learning, based on statistical emulation and Gaussian processes regression (Gaussian Process  Upper Confidence Bound optimisation, GPUCB
[43]).
5.1 Gaussian Processes  Upper Confidence Bound Optimisation
Gaussian Processes.
The key ingredient for the design problem is an efficient estimation of the unknown objective function, i.e. the average robustness as a function of the process kinetic parameters. Function approximation is a central task in machine learning and statistics. The general regression task can be formulated as follows
[8]: given a set of inputoutput pairs (training data), with and , determine a function s.t. is optimally close to the target values(usually in terms of minimising a suitable loss function). Several methods exist for addressing this task; in this paper we consider Gaussian Process (GP) regression, a popular Bayesian methodology
[41]. GPs are flexible nonparametric distributions over spaces of functions which can be used as prior distributions in a Bayesian framework, where the inputoutput pairs represent noisy observations of the unknown function. This enables a natural quantification of the uncertainty of the estimated function at every new input value; this uncertainty will play a central role in the optimal design strategy we propose in Section 5. We now give a semiformal definition of GP [41]:
Definition 3
A Gaussian Process over a (portion of) is a collection of random variables indexed by such that every finite dimensional marginal distribution is multivariate normal. Furthermore, there exist two functions (mean function) and (covariance function) such that the mean and covariance of the finite dimensional normal marginals is given by evaluating the mean and covariance functions at each point and each pair of points respectively.
We denote a sample from a GP with mean function and covariance function as
In practice, the inputoutput pairs in a regression task are often different features of experimentally observed data points. In this paper, the output points correspond to true functional evaluations of an unknown (and analytically intractable) function of the inputs. In this case, the regression task is often given the special name of emulation in the statistics literature: the true (but unknown) function is assumed to be a draw from a GP, and the functional evaluations are used as observations to obtain a posterior estimate of the unknown function. This approach was initially introduced in order to perform sensitivity analysis for deterministic computer models in [33]
; in that case, the function evaluations could be assumed to be noiseless (apart from numerical errors that were considered negligible in that paper). In our case, the function linking model parameters to average robustness cannot be computed, and we can only obtain a sampling approximation through a Statistical Model Checking procedure. This means that our function evaluations will be noisy; by virtue of the Central limit theorem we can assume that, provided sufficient samples were used for the SMC estimates, the noise in the observed robustness estimates will be approximately Gaussian.
^{1}^{1}1Note that here we approximate as a GP the average robustness score (or any other fitness score) as a function of parameters. We are not imposing any (Gaussian) approximation of the process itself. This therefore enables us to obtain an analytical estimate of the posterior process [41]. Furthermore, the SMC samples also allow us to estimate the (sample) variance in the average robustness at every sampled parameter value; this information can also be included leading us to a heteroscedastic (i.e. with non identical noise) regression problem (which is however still analytically tractable).
Gaussian Processes Optimization: the GPUCB algorithm
As we have seen, GP emulation provides a convenient way to explore approximately the average robustness of a stochastic process for different values of the model parameters. One could then be tempted to also use the emulated robustness profile for model design, i.e. find the optimum of the emulated function. This strategy, while appealing in its simplicity, is vulnerable to local optima: the emulated function is estimated based on relatively few function evaluation, so that, while the emulator typically provides a good approximation of the true function near the sampled points, regions of parameter space far from the sampled points may contain the true maximum undetected. Using the language of reinforcement learning, maximising the emulated function would privilege exploitation (i.e. using currently available information) at the expense of exploration. Obviously, given sufficient computational power, one may consider sampling many parameter points so as to have sufficient coverage of the whole region of interest; this strategy is however bound to fail in even moderate dimensions due to the curse of dimensionality.
An elegant solution to the above conundrum can be obtained by also considering the uncertainty of the emulated function (which is also computed analytically in GP regression): intuitively, one should explore regions where the maximum could plausibly be
, i.e. regions in parameter space where there is substantial posterior probability mass for the function to take a high value. We formalise these ideas in a recursive search rule, the so called Gaussian Process Upper Confidence Bound (GPUCB) algorithm: assume we have computed the average robustness at
parameter values (so that we have input output pairs). Let and be the mean and variance of the GP emulator at a given point in input space (recall that the marginal at any point will be Gaussian)^{2}^{2}2We now denote the input as to emphasize that they are the parameters of a stochastic process. We select the parameter value for the next function evaluation according to the following rule(5.3) 
where is a parameter. Thus, the next point for exploration does not maximise the emulated function, but an upper confidence bound at a certain confidence level specified by the parameter
(the quantile can be obtained by applying the inverse probit transform to the parameter).
[43] proved that this algorithm converges to the global maximum of the unknown function with high probability (which can be adjusted by varying the algorithm’s parameters).The primary difficulty in applying GPUCB is that, in order to be able to apply the rule, the emulated function must be computed at a large number of points; while this is obviously not as onerous as evaluating the true unknown function (as the emulator is known analytically), it may still be problematic for high dimensional parameter spaces. Nevertheless, the algorithm can be applied effectively for moderate sized parameter spaces (of the order of 10 parameters), and modular construction may be used to extend to higher dimensional systems [29, 6].
5.2 Experimental Results
Schlögl system.
We set up the experiment as follows.
We combine the robustness degree of the formula of Section 4.1 and the satisfaction probability in the systems design problem asking, at the same time, to maximize the robustness degree constraining the probability value to remain above 0.75.
We varied uniformly in , fixing all other parameters to the values of Table 1. We ran the GPUCB optimisation algorithm by first estimating the robustness degree and the satisfaction probability, using statistical model checking, for 30 points sampled randomly and uniformly from the parameter space, and then using the GPUCB strategy to estimate the maximum of the upper bound function in a grid of 200 points. If in this grid a point is found with a larger value than those of the observation points, we compute the robustness and satisfaction probability also for this new point, and add it to the observations (thus changing the GP approximation). Termination happens when no improvement can be made after three grid resamplings. Further integration of local maximisation can further improve the method.
In the experiment, repeated 10 times, we used a GP with radial basis kernel [8], with length scale fixed to 0.5 (after standardisation of the parameter range to
). The amplitude of the kernel was adaptively set to 60% of the difference between the maximum and the mean value of the robustness for the initial observations. The observation noise was experimentally fixed to 1, by monitoring the average standard deviation at different random parameter combinations.
Results are shown in table 2. As we can see, the result of the optimisation suggests that the more robust system satisfying the specification (i.e. remaining as much as possible above the threshold 300 for a sufficiently long amount of time) is the one obtained for . We can see that this is the case in Figure 5(b): the system becomes monostable, and stably remains above 550 units (corresponding to a robustness score above 250 with very high probability).
Parameter mean  Parameter range  Mean probability 

[979.31 999.99]  1  
Average Robustness  Number of function evaluations  Number of simulation runs 
348.97  34.4  3440 
Repressilator.
We consider a different optimisation problem, in which we keep model parameters constant and we try to optimise the parameters of the formula to make the robustness score as large as possible. This can be seen as a sort of dual problem, in which the goal is to learn the emergent behaviour of the model in terms of the most robustly satisfied formula (of fixed structure). Furthermore, the parametrisation of a formula is usually an underestimated problem, as the satisfaction/robustness heavily depends on these parameters. This problem has been partially tacked e.g. in [42] for deterministic models, but never for stochastic ones, to authors’ knowledge. In particular, we consider Formula (4.2) and optimise the temporal delays in the range and in the range . This can be seen as an attempt to learn the best bounds on the oscillatory period, through the filter of the logical specification of oscillations of equation (4.2).
In this experiment, we used the same settings of the optimisation algorithms as for the Schlóg system, save for the number of initial observations, set to 25 and constrained to lie in an equispaced grid. The parameters of the model are fixed to those shown in the caption of Figure 3(a).
In this case, due to the highly random duration of each oscillation cycle of the SHA model of the Repressilator, we expect to obtain a low value for and a large value for . This is indeed the case (see Table 3) confirming the intrinsic instability of the oscillatory pattern of the model. Moreover, the large variability in the value of can indicate a scarce contribution of the parameter in the determination of the robustness. We confirmed this intuition by emulating the average robustness score as a function of alone in the range for fixed to (data not shown), obtaining an essentially flat function: the average robustness varies between 15.01 (near ) and 6.89 (near ). We expect that an heteroschedastic treatment of noise could improve this estimate of .
Parameter mean  Parameter range  Mean probability 

[0, 500] [6963, 7000]  0.900  
Average Robustness  Number of function evaluations  Number of simulation runs 
17.20  35  3500 
6 Conclusion
Discussion.
In this paper we investigate a notion of robustness of behaviours of stochastic models, extending the robustness degree of STL formulae in a probabilistic setting. Discussing two case studies, a bistable model and the Repressilator, we showed that the distribution of the robustness degree of a formula provides more information than the satisfaction probability alone, and its average can be used to enforce robust behaviours by optimising it. Such optimisation is carried out using stateoftheart optimisation algorithms coming from reinforcement learning, which emulate the true function from just few samples, and perform very well in a simulation based scenario. Remarkably, the proposed approach to evaluate robustness and to system design can be applied both to CTMC and to SHA models. We also briefly considered the problem of learning the most effective parameters of a given formula, in the sense of finding the parameter combination maximising the robustness score. This hints towards a more ambitious goal, that of finding machine learning procedures to learn the emergent behaviours (described as temporal logic formulae) from models and from experimental data. Many problems need to be faced to achieve this goal, like how to learn formula structure, how to avoid overfitting (with respect to both formula structure and parameters), how to deal with the curse of dimensionality afflicting GPUCB and other optimisaton algorithms.
Related Work.
Temporal Logic (TL) is a very intuitive specification language to express formally the behavioural property emerging in a complex biological system. Several important extensions of TL, such as Metric Interval Temporal Logic (MITL) and Signal Temporal Logic (STL) [37], have been introduced to deal with densetime and realvalued signals, respectively.
In the last years there was a great scientific effort to enrich the classic qualitative semantics of TL or satisfiability (yes/no answer for the formula satisfaction of a trajectory) with more powerful and useful notions of quantitave semantics[27, 28, 24, 42, 3] (or robustness degree), providing a real value measuring the level of satisfaction or violation for a trajectory of the property of interest.
Several tools, such as BIOCHAM [16], STaLiRo [3] Breach [20], are now available to perform robustness analysis on the time series collected in wetlab experiments or produced by simulationbased techniques. The robustness degree have been successfully employed in the analysis of ODEbased biological models, to tune the parameters that discriminates the behaviours observed experimentally (the so called design problem).
Donze et al. in [22] proposed a multistep analysis, where they adopt STL to express dynamical properties and they use robustness and sensitivity analysis to sample efficiently the parameter space, searching for feasible regions in which the model exhibit a particular behavior. A similar method appeared also in a previous paper [21] of one of the coauthors.
In [6] the authors proposed a new approach, based on robustness degree, for the design of a synthetic biological circuit whose behaviour is specified in terms of signal temporal logic (STL) formulae. Also in this case stocasticity was not taken into account.
For what concern the stochastic models, while the satisfiability analysis has been considered as a discriminating criterion to tune the parameters in the design process using both simulationbased statistical approximated methods [14] and probabilistic exact methods [7], to the best of our knowledge, we are not aware of approaches using the robustness degree. Another related work in this sense is that of [36], where authors compute exactly upper and lower bounds on the satisfaction probability within a given region of the parameter space.
Future work.
The present work uses advanced machine learning concepts to address core problems in formal modelling; this is a relatively new line of work [14] which opens significant new avenues for further research. From the practical point of view, more extensive testing and an efficient and robust implementation (exploiting some of the possible parallelisms e.g. in SMC) will be important for the tool to be adopted. From the theoretical perspective, we plan to use multiobjective optimisation to find good parametrisation for conflicting objectives. Another interesting direction is to combine the design problem with the inference problem, which has recently been addressed for a number of continuous time stochastic systems [38]; this would open the possibility of addressing the control problem for such systems, simultaneously inferring the state of the system and designing the optimal input to lead it to a desired state.
Acknowledgements.
Work partially supported by the EUFET project QUANTICOL (nr. 600708) and by FRAUniTS.
References
 [1]
 [2] R. Alur, T. Feder & T.A. Henzinger (1996): The benefits of relaxing punctuality. J. ACM. Available at http://doi.acm.org/10.1145/227595.227602.
 [3] Y. Annapureddy, C. Liu, G. Fainekos & S. Sankaranarayanan (2011): STaLiRo: A Tool for Temporal Logic Falsification for Hybrid Systems. In: Proceedings of TACAS, doi:10.1007/9783642198359_21.
 [4] C. Baier, E.M. Clarke, V. HartonasGarmhausen, M.Z. Kwiatkowska & M. Ryan (1997): Symbolic Model Checking for Probabilistic Processes. In: Proc. of ICALP ’97, the 24th International Colloquium on Automata, Languages and Programming, Bologna, Italy, July 7 11, Lecture Notes in Computer Science 1256, Springer Berlin Heidelberg, pp. 430–440, doi:10.1007/3540631658_199.
 [5] C. Baier, B. Haverkort, H. Hermanns & J.P. Katoen (2003): ModelChecking Algorithms for ContinuousTime Markov Chains. IEEE Trans. Softw. Eng. 29(6), pp. 524–541, doi:10.1109/TSE.2003.1205180.
 [6] E. Bartocci, L. Bortolussi & L. Nenzi (2013): A temporal logic approach to modular design of synthetic biological circuits. In: In Proc. of CMSB 2013, the 11th International Conference on Computational Methods in Systems Biology, IST Austria, Klosterneuburg, Austria, September 2325, 2013, Lecture Notes in Computer Science 8130, SpringerVerlag, pp. 164–178, doi:10.1007/9783642391767.
 [7] E. Bartocci, R. Grosu, P. Katsaros, C. Ramakrishnan & S. A. Smolka (2011): Model Repair for Probabilistic Systems. In: Proceedings of TACAS 2011, the 17th International Conference on Tools and Algorithms for the Construction and Analysis of Systems, Lecture Notes in Computer Science 6605, Springer Berlin / Heidelberg, pp. 326–340, doi:10.1007/9783642198359_30.
 [8] C. M. Bishop (2006): Pattern Recognition and Machine Learning. Springer.
 [9] L. Bortolussi (2010): Limit behavior of the hybrid approximation of Stochastic Process Algebras. In: Proceedings of ASMTA 2010, doi:10.1007/9783642135682_26.
 [10] L. Bortolussi, J. Hillston, D. Latella & M. Massink (2013): Continuous Approximation of Collective Systems Behaviour: a Tutorial. Performance Evaluation 70(5), pp. 317–349, doi:10.1016/j.peva.2013.01.001.
 [11] L. Bortolussi & A. Policriti (2008): Hybrid approximation of stochastic process algebras for systems biology. In: Proceedings of IFAC WC, doi:10.3182/200807065KR1001.02132.
 [12] L. Bortolussi & A. Policriti (2010): Hybrid Dynamics of Stochastic Programs. Theoretical Computer Science 411(20), pp. 2052–2077, doi:10.1016/j.tcs.2010.02.008.
 [13] L. Bortolussi & A. Policriti (in print): (Hybrid) Automata and (Stochastic) Programs. The hybrid automata lattice of a stochastic program. Journal of Logic and Computation.
 [14] L. Bortolussi & G. Sanguinetti (2013): Learning and Designing Stochastic Processes from Logical Constraints. In: Proc. of QEST 2013, 10th International Conference on Quantitative Evaluation of Systems, Buenos Aires, Argentina, August 2730, 2013, 8054, pp. 89–105, doi:10.1007/9783642401961.
 [15] M. L. Bujorianu, J. Lygeros & M. C. Bujorianu (2005): Bisimulation for General Stochastic Hybrid Systems. In: Proceedings of HSCCl, pp. 198–214, doi:10.1007/9783540319542_13.
 [16] L. Calzone, F. Fages & S. Soliman (2006): BIOCHAM: an environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics 22, pp. 1805–1807, doi:10.1093/bioinformatics/btl172.
 [17] T. Chen, M. Diciolla, M. Kwiatkowska & A. Mereacre (2011): Timebounded verification of CTMCs against realtime specifications. In: Proc. of FORMATS 2011, the 9th International Conference on Formal Modeling and Analysis of Timed Systems, Aalborg, Denmark, September 21–23, Lecture Notes in Computer Science 6919, Berlin, Heidelberg, pp. 26–42, doi:10.1007/9783642243103_4.
 [18] M.H.A. Davis (1993): Markov Models and Optimization. Chapman & Hall.
 [19] Andrea Degasperi & Stephen Gilmore (2008): Sensitivity analysis of stochastic models of bistable biochemical reactions. In: Formal Methods for Computational Systems Biology, Lecture Notes in Computer Science 5016, Springer, pp. 1–20, doi:10.1007/9783540688945_1.
 [20] A. Donzé (2010): Breach, A Toolbox for Verification and Parameter Synthesis of Hybrid Systems. In: Proceedings of CAV. Available at http://dx.doi.org/10.1007/9783642142956_17.
 [21] A. Donzé, G. Clermont & C.J. Langmead (2010): Parameter Synthesis in Nonlinear Dynamical Systems: Application to Systems Biology. Journal of Computational Biology 17(3), pp. 325–336, doi:10.1007/9783642020087_11.
 [22] A. Donzé, E. Fanchon, L. M. Gattepaille, O. Maler & P. Tracqui (2011): Robustness analysis and behavior discrimination in enzymatic reaction networks. PLoS One 6(9), p. e24246, doi:10.1371/journal.pone.0024246.
 [23] A. Donzé, T. Ferrer & O. Maler (2013): Efficient Robust Monitoring for STL. In: Proc. of CAV 2013, the 25th International Conference on Computer Aided Verification, Saint Petersburg, Russia, July 1319, Lecture Notes in Computer Science 8044, pp. 264–279, doi:10.1007/9783642397998_19.
 [24] A. Donzé & O. Maler (2010): Robust satisfaction of temporal logic over realvalued signals. In: Proc. of FORMATS 2010, the 8th International Conference on Formal Modeling and Analysis of Timed Systems, Klosterneuburg, Austria, September 8–10, 6246, pp. 92–106, doi:10.1007/9783642152979_9.
 [25] R. Durett (2012): Essentials of stochastic processes. Springer, doi:10.1007/9781461436157.
 [26] M.B. Elowitz & S. Leibler (2000): A synthetic oscillatory network of transcriptional regulators. Nature 403, pp. 335–338, doi:10.1038/35002125.
 [27] G. Fainekos & G. Pappas (2007): Robust Sampling for MITL Specifications. In: Proc. of FORMATS 2007, the 5th International Conference on Formal Modeling and Analysis of Timed Systems, Lecture Notes in Computer Science 8044, pp. 264–279, doi:10.1007/9783540754541_12.
 [28] G. E. Fainekos & G. J. Pappas (2009): Robustness of temporal logic specifications for continuoustime signals. Theor. Comput. Sci. 410(42), pp. 4262–4291, doi:10.1016/j.tcs.2009.06.021.
 [29] A. Georgoulas, A. Clark, A. Ocone, S. Gilmore & G. Sanguinetti (2012): A subsystems approach for parameter estimation of ODE models of hybrid systems. In: Proc. of HSB 2012, the 1st International Workshop on Hybrid Systems and Biology, EPTCS 92, pp. 30–41, doi:10.4204/EPTCS.92.3.
 [30] D.T. Gillespie (1977): Exact Stochastic Simulation of Coupled Chemical Reactions. J. of Physical Chemistry 81(25), doi:10.1021/j100540a008.
 [31] R. Gunawan, Y. Cao, L. Petzold & F.J. Doyle III (2005): Sensitivity analysis of discrete stochastic systems. Biophysical Journal 88(4), p. 2530, doi:10.1529/biophysj.104.053405.
 [32] K. D. Jones, Konrad V & D. Nickovic (2010): Analog property checkers: a DDR2 case study. Formal Methods in System Design 36(2), pp. 114–130, doi:10.1007/s107030090085x.
 [33] M. Kennedy & A. O’Hagan (2001): Bayesian Calibration of Computer Models. Journal of the Royal Stat. Soc. Ser. B 63(3), pp. 425–464, doi:10.1111/14679868.00294.
 [34] M. Komorowski, M. J. Costa, D. A. Rand & M. PH Stumpf (2011): Sensitivity, robustness, and identifiability in stochastic chemical kinetics models. PNAS USA 108(21), pp. 8645–8650, doi:10.1073/pnas.1015814108.
 [35] Marta Kwiatkowska, Gethin Norman & David Parker (2004): Probabilistic symbolic model checking with PRISM: a hybrid approach. Int. J. Softw. Tools Technol. Transf. 6(2), pp. 128–142, doi:10.1007/s1000900401402.
 [36] S. Drazan L. Brim, M. Ceska & D. Šafránek (2013): Exploring Parameter Space of Stochastic Biochemical Systems using Quantitative Model Checking. In: Proc. of CAV 2013, the 25th International Conference on Computer Aided Verification, Saint Petersburg, Russia, July 1319, Lecture Notes in Computer Science 8044, pp. 107–123, doi:10.1007/9783642397998_7.
 [37] O. Maler & D. Nickovic (2004): Monitoring Temporal Properties of Continuous Signals. In: Proc. of Joint International Conferences on Formal Modeling and Analysis of Timed Systmes, FORMATS 2004, and Formal Techniques in RealTime and Fault Tolerant Systems, FTRTFT 2004, Grenoble, France, September 2224, 3253, pp. 152–166, doi:10.1007/9783540302063_12.
 [38] A. Ocone, A. J. Millar & G. Sanguinetti (2013): Hybrid regulatory models: a statistically tractable approach to model regulatory network dynamics. Bioinformatics 29(7), pp. 910–916, doi:10.1093/bioinformatics/btt06.
 [39] M. Opper, A. Ruttor & G. Sanguinetti (2010): Approximate inference in continuous time GaussianJump processes. In: Proceedings of NIPS 2010, the 4th Annual Conference on Neural Information Processing Systems, 6–9 December 2010, Vancouver, British Columbia, Canada, pp. 1831–1839. Available at http://books.nips.cc/papers/files/nips23/NIPS2010_1095.pdf.
 [40] Amir Pnueli (1977): The temporal logic of programs. Foundations of Computer Science, IEEE Annual Symposium on 0, pp. 46–57, doi:10.1109/SFCS.1977.32.
 [41] C. E. Rasmussen & C. K. I. Williams (2006): Gaussian Processes for Machine Learning. MIT Press.
 [42] A. Rizk, G. Batt, F. Fages & S. Soliman (2008): On a Continuous Degree of Satisfaction of Temporal Logic Formulae with Applications to Systems Biology. In: Proc. of CMSB 2008, the 6th International Conference on Computational Methods in Systems Biology, Rostock, Germany, October 12–15, Lecture Notes in Computer Science 5307, pp. 251–268, doi:10.1007/9783540885627_19.
 [43] Niranjan Srinivas, Andreas Krause, Sham M. Kakade & Matthias W. Seeger (2012): InformationTheoretic Regret Bounds for Gaussian Process Optimization in the Bandit Setting. IEEE Transactions on Information Theory 58(5), pp. 3250–3265, doi:10.1109/TIT.2011.2182033.
 [44] H. L. S. Younes, M. Z. Kwiatkowska, G. Norman & D. Parker (2004): Numerical vs. Statistical Probabilistic Model Checking: An Empirical Study. In: Proc. of 2004, the 10th International Conference on Tools and Algorithms for the Construction and Analysis of Systems, Barcelona, Spain, March 29  April 2, doi:10.1007/9783540247302_4.