Parameter Estimation of Biological Pathways Using Data Assimilation and Model Checking

Abstract

This paper presents a novel method to estimate kinetic parameter of biological pathways by using observed time-series data and other knowledge that cannot be formulated in the form of time-series data. Our method utilizes data assimilation (DA) framework and model checking (MC) technique, with a quantitative modeling and simulation architecture named hybrid functional Petri net with extension (HFPNe). Proposed method is applied to an HFPNe model underlying circadian rhythm in mouse. We first translate 23 rules of biological knowledge with temporal logic for the model checking, which are not described in the time-series data. Next, we employ particle filter often applied to DA for our estimation procedure. Each particle checks whether its simulation result satisfies the rules or not, and the result of the checking is used for its resampling step. Our simulation results show that proposed method is faster and more accurate than previous method.

Keywords

• Hybrid functional Petri net with extension
• parameter estimation
• data assimilation
• particle filter
• model checking
• temporal logic

-

Modeling and simulating large-scale biological pathways have played an important role in systems biology. Owing to their importance, many formal description methods of biological pathway models have been made so far [ 1–3 ]. Petri net and its related concepts are one of the succeeded ways of describing biological models [ 4–6 ], which have been used for modeling a wide variety of biological pathways and succeeded in reproducing consistent time-series profiles of biological elements such as the concentrations of mRNAs and proteins by means of computer simulations.

Simulation studies on biological pathways promises a deep understanding of complex cellular mechanisms by investigating the dynamic feature. Simulationbased models are commonly governed by a series of parameters, e.g. initial values, reaction speeds and threshold values of activities. Before the model can be simulated, all parameters must be assigned in advance. However, most parameters are often unknown or not obvious. In general, such parameters are carefully tuned by experts to fit the simulated elements with observed in vivo/vitro experiment results. Due to the nonlinearity of the model, parameter estimation is difficult and requires a lot of trial and errors. Small differences of the parameters make large gap between reality and simulation results. Therefore, conventional hand tuning method severely limits the size and complexity of simulation models built as more output data (e.g., microarray gene expression data) are being measured.

The aim of this paper is to develop a novel method to automatically and efficiently estimate kinetic parameters of a given model or a model starting from scratch by combining data assimilation (DA) and model checking (MC) approaches, coupled with observed experiment data. Observed experiment data includes well-defined time-series data and other knowledge that cannot be formulated in the form of time-series data.

DA was originally established in the field of geophysical simulation science. Nagasaki et al. [ 7 ] have proposed a so-called genomic data assimilation approach. Their DA framework enables users to handle both the model construction and parameter tuning in the context of statistical inferences, and establishes a link between the HFPNe simulation model and observed data, e.g., microarray gene expression data [ 7 ] or time series proteomic data [ 8 ]. Current DA approach has some issues because it depends on providing successive time points (at least 1020 time points) of time-series data by biological experiments. That is, for a small time-series data set including for example two or three time points, DA will cost massive computational resource, in some cases, it will be completely impossible to estimate parameters. The response to this difficulty of dealing with sparse and/or not well-defined time-series data is the use of model checking.

MC is a method for automatic verification of system requirements [ 9 ], which firstly used in hardware field because of its determined behavior and limited value space. Today, this technique has been applied to more complex biological models. There are several studies that use model checking for parameter estimation, e.g. Donaldson and Gilbert developed a computational system named MC2(GA) [ 10 ] that couples model checking with genetic algorithm for parameter estimation. Li et al. [ 11 ] proposed online model checking approach based parameter estimation framework applied to the HFPNe class. In order to check the model, one need to write biological rules of the knowledge with temporal logic. For example, “A biological phenomena F oo keeps decreasing until Bar rises.” can be written as “d([Foo])≤0 U d([Bar])>0” in a kind of temporal logic. Various knowledge can be described in this way. Nevertheless, in order to improve the estimation efficiency and accuracy, it is expected to find a general methodology to determine parameters by combining DA and MC dealing with both time-series experimental data and biological queries.

The paper is organized as follows. In Methods, we briefly explain how to (i) construct biological models with HFPNe, (ii) estimate unknown parameters in a nonlinear state space model, (iii) translate biological rules with a temporal logic for querying system properties, and (iv) combine MC with DA for parameter estimation. In Results, we compare our novel method with previous method by applying them to mouse circadian rhythm model represented by HFPNe. We show that our method is potentially faster and more accurate than previous one that excludes MC technique. 2 2.1

Methods

Hybrid Functional Petri Net with Extension (HFPNe)

HFPNe is developed as a biosimulation tool for pathway modeling and simulation extended from original Petri net [ 13 ]. HFPNe can deal with three types of data: discrete, continuous and generic, whereas the original Petri net deal with only discrete data. HFPNe consists of three types of elements: entity, process and connector (see Figure 1 (a)). Figure 1 (b) shows connection rules in HFPNe. For more definitions and usages of HFPNe, see Nagasaki et al. [ 13 ].

Figure 2 shows circadian rhythm model of mouse represented by HFPNe [ 14 ]. This model is composed of of 12 entities, 28 processes and 45 connectors. Due to the space restriction, Appendix gives the details of these elements. Initial value of entities mi(0) (i=1, · · ·, 12), reaction speed parameters ki (i=1, 2 and ki is a common parameter to control speeds of similar biological processes: k1 for protein binding; k2 for translation), and threshold parameters si (i=1, · · ·, 3) are unknown parameters. 2.2

Data Assimilation for Parameter Estimation

We here explain how to estimate parameters in a simulation model from timeseries data with the use of DA.

Data Assimilation with Nonlinear State Space Model DA is an approach to improve the accuracy of the models by combining with data, which can deal with models formulated by nonlinear state space model (SSM) given by two equations [ 15 ]: mt=f (mt−1, wt, θsys) , (1) yt=Hmt + t . (2) Equation (1) is called “system model” and Equation (2) is called “observed model”. In the system model, mt≡(m1t, · · ·, mpt)T is a state vector consisting of p state variables mit (i=1, · · ·, p) at discrete time point t. wt denoting system noise is an l-dimensional white noise at t with a density q(w). f is a vector-valued function, f : Rp+l → Rp, and θsys is a vector of model parameters. In the observed model, yt∈Rd is an observation vector at t, H∈Rd × Rp is an observation matrix; Hij takes value one if observed value of jth entity corresponds to the ith element of yt, otherwise zero. t called observational noise is a d-dimensional white noise at t with a density r( ). The likelihood of the parameter is given by

L(θsys)=p(y1, · · ·, yT |θsys)= tT=1 p(yt|Y t−1, θsys)= tT=1 p(yt|mt, θsys)p(mt| Y t−1, θsys)dmt, where Y t={y1, · · ·, yt}. The maximum likelihood estimator (MLE) for θsys is given by argmax log L(θsys).

θsys We can estimate parameters with MLE, however, this maximum likelihood method has an important issue that the value of log L(θsys) computed by Monte Carlo filter includes an approximation error. For accurate estimation, huge computation resources are thus required.

Self-Organizing State Space Model To deal with the difficulty of maximum

likelihood method, we use self-organizing state space model (SOSSM) [ 16–18 ]. We can estimate parameters based on Bayesian inference with SOSSM by weaving parameters in the state vector as The SSM for this vector is given by zt=F ∗(zt−1, wt) and yt=H∗zt+ t, where zt= mt θsys

.

F ∗(zt−1, wt) = f (mt−1, wt, θsys)

θsys H∗zt = Hmt. We can obtain marginal posterior densities without obtaining MLE of θ as p(mT | Y T ) =

p(zT | Y T )dθsys p(θsys | Y T ) =

p(zT | Y T )dmT .

In this way, we can avoid the issue of maximum likelihood method during the estimation.

Particle Filter As mentioned above, we have to calculate distribution p(zT |Y T ) for estimating parameters. However, it generally becomes a non-gaussian distribution in SSM. Therefore, it is needed to represent this distribution with some method. We here use a sequential Monte Carlo method called particle filter (PF) [ 19 ]. Figure 3 shows the overview of particle filter’s algorithm. In particle filter, predictive distribution p(zt|Y t−1) and filter distribution p(zt|Y t) are approximated by m in which each realization is called particle as follows: pt|t−1 ≡ {pt(|1t)−1, · · · , pt(|mt−) 1} ∼ p(zt | Y t−1)

pt|t ≡ {pt(|1t), . . . , pt(|mt )} ∼ p(zt | Y t), where pt(|jt)−1 (j=1, · · ·, m) and pt(jt) (j=1, · · ·, m) are (p+|θsys|)-dimensional num| bers. This algorithm is processed with the following steps. Step 1 Generate the (p+|θsys|)-dimensional random number p(0j0) for j=1, · · ·, m. | Step 2 Repeat the following three steps for the observed time points t=1, · · ·, T .

Step 2.1 (Prediction step) Generate m system noises wt(j) independently and identically from q(wt), and compute particles by inputting filtered states into simulation model, pt(|jt)−1=F ∗(pt(−j)1|t−1, wt(j)) for j=1, · · ·, m.

Step 2.2 (Weight calculation step) Compute the weights of importance

for the particles according to

αt(j)=p(yt|pt(|jt)−1)=r(yt − H∗pt(|jt)−1) for j=1, · · ·, m. These weights are then normalized as α¯t(j)=αt(j) Step 2.3 (Resampling step) Generate filtered state pt(jt) for j=1, · · ·, m | by resampling {pt(|1t)−1, · · · , pt(|mt−) 1} with the probabilities {α¯t(1), · · · , α¯t(m)}. (3) m j=1αt(j). 2.3

Model Checking

We explain how to represent requirements of biological pathway models for model checking. Model checking is a method to verify whether models satisfy the requirements or not. Temporal logic formulae are often used to describe the system requirements. In this study, we selected PLTL (Probabilistic Linear-time Temporal Logic) for querying dynamic models of cellular networks [ 20, 21 ] which extends original LTL to a stochastic setting with a probability operator and a filter criterion defining the starting state where the property is satisfied. [PLTL Syntax] Table 1 shows definition of PLTL syntax, which is used to ask for the probability of user’s query via a PLTL formulae ψ. In the LTL expression φ{AP }, φ will be checked from the state that AP is satisfied rather than from the default initial state, where AP is called atomic proposition and takes boolean domain. PLTL allows (i) LTL expression to contain temporal operators, i.e., X, F, G, U, R. Five temporal operators are used to describe the sequencing of the states along the execution; and (ii) the usage of ψ without probabilistic operators (i.e. simply in the form of LTL), which is useful when the model is deterministic. [PLTL Semantics] The semantics of PLTL is defined over the finite sets of finite paths through system’s state space, obtained by repeated simulation runs of HFPNe models. The PLTL formula is built upon two components: probabilistic operator and property LTL. For each simulation run, the LTL expression is evaluated to a boolean truth value, and the probability of the LTL statement holding true is calculated based on the whole set of simulation results.

For the probability operator components, there are two distinct operators: (i) P x(LT L) is any inequality comparison of the probability of the property LTL holding true, for example P≥0.5(LT L); and (ii) P=?(LT L) returns the value of the probability of the property holding true. The semantics of the temporal logic operators are described in Table 2. Concentrations of biochemical species in the model are denoted by [variableName]. A special variable, [time], stands for simulation time.

Due to the ability of PLTL, it is possible to define functions of two different natures: functions that return a real number and functions that return a boolean value. An example of the real number function is d([variableN ame]) which returns the subtracted value of [variableName] between time i and i−1. Note that, d([variableN ame]) equals zero at time point zero. One example of a boolean function is similarAbsolute(value a, value b, value ), which returns true if |a−b|< or else it returns false. Table 3 shows the rules written in PLTL for circadian rhythm model of mouse. As stated in Introduction, we have applied either DA or MC to pathway model in the previous researches, but used time scales are different. We here employ common time scale called Petri net time for combining DA with MC, which is the virtual time unit of the HFPNe model denoted by [pt]. We define: (i) simInt∈R as simulation interval; (ii) mcInt∈R which is a multiple of simInt as a model checking interval; and (iii) M apOttoP t : N→R as a mapping from observed time point to Petri net time for combining. From the MC’s viewpoint, Xφ means that φ must be true at the state after mcInt Petri net time. Meanwhile, a simulation from time M apOttoP t(t−1) to M apOttoP t(t) is a run in f from DA’s viewpoint.

Hu¨rseler and Ku¨nsch mentioned that it is difficult to generate a good initial distribution of parameters with SOSSM [ 22 ]. This is because parameters in resampled particles are the subset of parameters in the initial particles and a model with randomly generated parameters rarely satisfy all rules. To generate good initial distribution of parameters, two considerations are designed:

Rule LTL translation

Rule 1 GC(o[npceernmtrRaNtAi]on< o3f.p8er&&m[RpNerAmRisNAb]et>we0e.n2)0.2 and 3.8.

Rule 2 Concentration of Rev-Erv mRNA is between 0.2 and 3.4.

G([Rev-Erv mRNA] < 3.4 && [Rev-Erv mRNA] > 0.2) Rule 3 Concentration of Bmal mRNA is between 0.2 and 4.3.

G([Bmal mRNA] < 4.3 && [Bmal mRNA] > 0.2) Rule 4 Concentration of Bmal mRNA is between 2.4 and 2.6 after time becomes 20.

G(similarAbsolute([Clock mRNA],2.5,0,1)){[time] > 20} Rule 5 Concentration of Cry mRNA is between 3.8 and 0.2.

G([Cry mRNA] < 3.8 && [Cry mRNA] > 0.2) Rule 6 GC(o[nPcEeRn]tr<at2i.o6n &o&f P[EPERR]is>b0e.tw2)een 2.6 and 0.2.

Concentration of CRY is between 2.6 and 0.2.

Rule 7 G([CRY] < 2.2 && [CRY] > 0.2)

Concentration of PER/CRY is between 2.7 and 0.2.

Rule 8 G([PER/CRY] < 2.7 && [PER/CRY] > 0.2) Rule 9 Concentration of REV ERB is between 1.5 and 0.2.

G([REV ERB] < 1.5 && [REV ERB] > 0.2) Rule 10 GL(odc(a[lpmeraxmRimNAu]m) ≥con0ce&n&trXa(tdio(n[poefr pmeRrNAm])R<N0)A ⇒is g[rpeeartemrRNtAh]an> 22.0..0) Rule 11 GL(odc(a[lpmerinmiRmNuA]m)≤co0n&c&enXt(radt(i[opnerofmRpNeAr]m)>R0)N⇒A [ispelerssmRtNhAa]n <1.10..0) Rule 12 GL(odc(a[lRmeva-xEirmvummRNAc]o)n≥ce0nt&r&atXio(nd(o[fReRve-vE-rEvrmvRNmAR])N<0A) is⇒greater than 1.5. [Rev-Erv mRNA] > 1.5) Rule 13 GL(odc(a[lRmevi-nEimrvummRNcAo]n)c≤e0nt&ra&tiXo(nd(o[fRReve-vE-rEvrvmRmNAR])N>A0)is less than 1.0. ⇒ [Rev-Erv mRNA] < 1.0) Rule 14 GL(odc(a[lBmmaalxmimRNuAm])c≥on0ce&n&trXa(tdio(n[BomfalBmmRaNlA]m)<R0N)A⇒is[gBrmeaaltemrRNthA]an>11.5..5) Rule 15 GL(odc(a[lBmmailnimmRuNAm])c≤on0c&e&ntXra(tdi(o[nBmoaflBmmRNaAl]m)>R0)NA⇒is[BlemsasltmhRaNnA]1.<0.1.0) Rule 16 GL(odc(a[lCmryaxmRimNAu]m)≥c0on&c&enXt(rda(t[ioCnryomfRCNAry])m<0R)N⇒A i[sCrgyremaRtNeAr]th>an2.20.)0.

Rule 17 GL(odc(a[lCmryinmiRmNuA]m)≤co0n&c&enXt(radt(i[oCnryofmRCNrAy])m>0R)N⇒A i[sClreyssmtRhNaA]n 1<.01..0) Rule 18 tWhahnencocnocnecnetnrtartaitoinonofopfeBrmmaRlNmAR.NA takes local minimum, concentration of Bmal mRNA is less

G(d([Bmal mRNA])≤0 && X(d([Bmal mRNA])>0) ⇒ [Bmal mRNA] < [per mRNA]) Rule 19 When concentration of Bmal mRNA takes local maximum, concentration of Bmal mRNA is grater than concentration of per mRNA

G(d([Bmal mRNA])≥0 && X(d([Bmal mRNA])<0) ⇒ [Bmal mRNA] > [per mRNA]) Rule 20 cWonhceenntcroantcieonntroaftpioenr mofRpNerAm.RNA takes local minimum, concentration of per mRNA is less than

G(d([per mRNA])≤0 && X(d([per mRNA])>0) ⇒ [per mRNA] < [Bmal mRNA]) Rule 21 tWhahnencocnocnecnetnrtartaitoinonofoBfmpearl mmRRNNAA takes local maximum, concentration of per mRNA is grater

G(d([per mRNA])≥0 && X(d([per mRNA])<0) ⇒ [per mRNA] > [Bmal mRNA]) Rule 22 GC(o[npceernmtrRaNtAi]on> o[fRpeevr-EmrvRmNRANA]is) greater than concentration of Rev-Erv mRNA.

Rule 23 GC(o[nCcLeOnCtKr]at>io[nPEoRf]C)LOCK is greater than concentration of PER.

Firstly, we repeat particle filters many times by regenerating p0|0 from pT |T like a crossover in genetic algorithm since the length of time courses for biological use is generally shorter than that for other fields; Secondly, a threshold value T h is used to discard particle whose unsatisfied number by the checking is greater than T h. T h is changed for each run of particle filter. In this study, system noise is not taken into account in order to accelerate the estimation. Nevertheless, a generic process can be mapped to Java object supporting HFPNe model in a nondeterministic settings.

Algorithm 1 shows pseudocode of our parameter estimation method. In step 9, P redictandM C returns predicated particle and the result whether it is worthless or not. The detail of this function is displayed in Algorithm 2. In steps 10–15, if the simulation results unsatisfy more than T h rules, the weight of the particle will become zero, or else the weight is calculated by Equation (3) via r. To calculate r, we assume that observed noise is a gaussian white noise and its mean is zero. The variance of observed noise is thus needed to be estimated. We generate multiple candidate values and use the value that has the maximum likelihood as the variance. Steps 17–19 are designed to break the run of particle filter if all the particles unsatisfied more than T h rules; Otherwise T h is incremented or decremented in steps 23–27. p0|0 for the next run of particle filter is generated from pT |T in step 28.

Algorithm 1. Pseudecode of our parameter estimation method.

In Algorithm 2, steps 2–3 convert time scale from observed time to Petri net time. Step 4 separates particle p into m(start) and θ, where m(t) is a vector consisting of the values of p entities at Petri net time t and θ is a vector of parameters to be estimated. m(pt+simInt) is predicted by simulation in step 6. Step 7 returns the number of unsatisfied rules via M odelChecking(). In steps 9–11, if simulation result violates more than T h rules, simulation will stop immediately.

Algorithm 2 Prediction step of particle filter and model checking. 1: function P redictandMC(p, t, Rules, T h) 2: start ← MapOttoP t(t − 1) 3: end ← MapOttoP t(t) − simInt 4: (m(start), θ) ← p 5: for pt ← start to end step simInt do 6: m(pt + simInt) ← Simulation(m(pt), θ) 7: if (pt + simInt)%mcInt = 0 then 8: N ← ModelChecking(m, Rules) 9: if T h < N then 10: return ([m(pt + simInt)T , θT ]T , false) 11: end if 12: end if 13: end for 14: return ([m(MapOttoP t(t))T , θT ]T , true) 3

Results and Discussions

3.1

Estimation environment and evaluation criteria

We estimate 17 parameters of circadian rhythm model of mouse (i.e., 12 initial values of entities, two reaction speeds and three threshold values). We use synthesized data set coupled with simulation and observed noise t ( t ∼ N (0, 0.052I12)) as observed data. It contains 312 data of 26 time points (see Figure 4) for each 12 biological entities – five mRNAs, five proteins and two complex proteins –. One observed time point is mapped to five Petri net times. Table 4 summarizes details of our estimation environment. Parameter search range is set from zero to 15.

To evaluate the results of estimation, following Scores are defined and calculated before the step RegenerateP articles step given in Algorithm 1.

Scoremean = tT=1 ep=1 |SimResult(P ar|YamTsm|ean,e,t)−yt,p|2 Scoremode = tT=1 ep=1 |SimResult(P ar| YamTs|mode,e,t)−yt,p|2 Scoremedian = tT=1 ep=1 |SimResult(P ara|YmsTm|edian,e,t)−yt,p|2 Scorebest = tT=1 ep=1 |SimResult(P a|rYamTs|best,e,t)−yt,p|2

Score = min{Scoremean, Scoremode, Scoremedian, Scorebest, Scorecurrent}, where yt,p is an observed value of entity p at observed time point t. P aramsbest denotes parameters of a particle which is made of particles. P aramsmean, P aramsmode and P aramsmedian represent mean, mode and median value of all particles’ parameters respectively. SimResult(P arams, e, t) returns the value of entity e at time t which is calculated by simulation with P aram. Scorecurrent is ∞ at first time of the calculation, otherwise it is the Score of previous calculation. 3.2

Experimental Result

Comparison between PFMC and PF We compare the performance between our new method Particle Filter with Model Checking (PFMC) and previous method Particle Filter (PF). The estimation experiments are carried out on workstation of Intel Xeon E5450 (3.0GHz) with 32G bytes of memory. We performed estimation 100 times for both methods. Figure 5 shows the result with respect to the distribution of Score with elapsed time. Mean of Score is also analyzed by Welch’s t-test (See Table 5). Nearly all Score of PFMC and PF are good on and after 600 seconds. Both medians are also good on and after 600 seconds, but there are many bad cases before 1,000 seconds for PF. That is, roughly speaking, if there are enough amounts of observed data, there is no much difference by using either PMFC or PF for the estimation. However, in almost all cases, we can finish the estimation within a short time incorporating model checking.

Estimation with small amount of observed data To investigate the per

formance in the case of small amount of observed data, we use only first ten time points of five biological entities’ data for the estimation: per mRN A, Cry mRN A, Rev Erv mRN A, Clock mRN A and Bmal mRN A. More detailedly, we use all default parameters for Score calculation and just overwrite p to five and T to ten. The estimation results are exhibited in Figure 6 and Table 6.

The results clearly show that on and after 1,200 seconds, in contrast to the previous experiment, there is difference not only between bad cases, but also between medians of PFMC and PF. This is because it is difficult to estimate parameters which makes certain rhythms with only two cycles of observed data. Therefore, median Scores of PF are not good before 3,000 seconds. Moreover, convergence of PFMC is worse than previous experiment.

Effects of the rules We also investigate the effect of rules by checking which rule is unsatisfied which results in the cutting of a particle. We run estimation 100 times with default parameters and all of observed data. All the runs finished after 20 minutes and the results are shown in Figure 7.

Two kinds of effects used in model checking approach can be considered. First, it cuts useless particles,which enables us to estimate with more particles or more runs of particle filters. Second, it cuts bad results which facilitates the estimation only with observed data. From the first effect’s viewpoint, the rule that cut more particles is a good rule. This is due to the fact that rule is able to cut many particles from early time because the number of unique particle decreases with the time elapse in our method. Rule 1 to 3, 5 to 9, 22 and 23 are such rules. From the second effect’s viewpoints, good rules are different depending on behaviors of target models. For the circadian rhythm model, it is important to reproduce the oscillations within a certain range. Rule 10 to 21 are specified for verifying the behavior of oscillation. Generally, it is not easy to prepare this kind of rules before trial. Nevertheless, unlike observed data, it will be a great help in improving the efficiency and accuracy of conventional parameter estimation process and eventually leading to better understanding of biological pathways. 4

Conclusions

We propose a novel parameter estimation method for biological pathways. By combining model checking with DA framework, our method enables us to use various knowledge in addition to observed time series data. We extract 23 biological rules with temporal logic which cannot be formulated in the form of time-series data. Proposed method and previous method are applied to mouse circadian rhythm model of HFPNe by means of performance evaluation. Results shows that (i) if estimations execute with enough amounts of observed data, our method can practically give good parameters in a short time; and (ii) if estimations execute with small amount of observed data, new method is much faster than the method without model checking.

Appendix:

Corresponding connector c45 c40 c39

References

1. Yan , H. , Zhang, B. , Li , S. , Zhao , Q. : A Formal Model for Analyzing Drug Combination Effects and its Application in TNF-alpha-induced NFkappaB Pathway . BMC Syst Biol ., 4 ( 50 ) ( 2010 )

2. Antoniotti , M. , Policriti , A. , Ugel , N. , Mishra , B. : Model Building and Model Checking for Biochemical Processes . Cell Biochem . Biophys., 38 ( 3 ), 271 - 286 ( 2003 )

3. Regev , A. , Silverman , W. , Shapiro , E. : Representation and Simulation of Biochemical Processes Using the Pi-calculus Process Algebra . Pac. Symp. Biocomput., 459 - 470 ( 2001 )

4. Chaouiya , C. : Petri Net Modelling of Biological Networks . Brief Bioinform., 8 ( 4 ), 210 - 219 ( 2007 )

5. Koch , I. , Heiner , M. : Petri nets . In Analysis of Biological Networks . Edited by Junker, B.H. , Schreiber , F. A Wiley Interscience Publication, 139 - 180 ( 2008 ).

6. Matsuno , H. , Li , C. , Miyano , S. : Petri Net Based Descriptions for Systematic Understanding of Biological Pathways . IEICE Trans. Fundamentals, E89-A(11) , 3166 - 3174 ( 2006 )

7. Nagasaki , M. , Yamaguchi , R. , Yoshida , R. , Seiya , I. , Doi , A. , Tamada , Y. , Matsuno , H. , Miyano , S. : Genomic Data Assimilation for Estimating Hybrid Functional Petri Net from Time-course Gene Expression Data . Genome Inform., 17 ( 1 ), 46 - 61 , 2006 .

8. Tasaki , S. , Nagasaki , M. , Kozuka-Hata , H. , Semba , K. , Gotoh , N. , Hattori , S. , Inoue , J. , Yamamoto , T. , Miyano , S. , Sugano , S. , Oyama , M. : Phosphoproteomics-based Modeling Defines the Regulatory Mechanism Underlying Aberrant EGFR Signaling . PLoS One , 5 ( 11 ), e13926 ( 2010 )

9. Be´rard, B. , Bidoit , M. , Finkel , A. , Laroussinie , F. , Petit , A. , Petrucci , L. , Schnoebelen , P. , McKenzie , P. : Systems and Software Verification: Model-Checking Techniques and Tools , Springer ( 2001 )

10. Donaldson , R. , Gilbert , D. : A Model Checking Approach to the Parameter Estimation of Biochemical Pathways . In: Proceedings of the 6th International Conference on Computational Methods in Systems Biology (CMSB '08) , pp. 269 - 287 . SpringerVerlag, Berlin, Heidelberg ( 2008 ).

11. Li , C. , Nagasaki , M. , Hock Koh , C. , Miyano , S. : Online model checking approach based parameter estimation to a neuronal fate decision simulation model in C. elegans with hybrid functional Petri net with extension , Mol. Biosyst ., 7 ( 5 ), 1576 - 1592 ( 2011 )

12. Nagasaki , M. , Saito , A. , Jeong , E. , Li , C. , Kojima , K. , Ikeda , E. , Miyano , S. : Cell Illustrator 4 .0: A Computational Platform for Systems Biology , In Silico Biol., 10 , 0002 ( 2010 )

13. Nagasaki , M. , Doi , A. , Matsuno , H. , Miyano , S.: A Versatile Petri Net Based Architecture for Modeling and Simulation of Complex Biological Processes . Genome Inform., 15 ( 1 ), 180 - 197 ( 2004 )

14. Circadian rhythms in Mus musculus , http://www.csml.org/models/ csml-models/ circadian-rhythms-in-mouse/

15. Kitagawa , G.: Non-Gaussian State-space Modeling of Nonstationary Time Series . Journal of the American Statistical Association , 82 ( 400 ), 1032 - 1063 , 1987 .

16. Kitagawa , G.: Self-organizing State Space Model . J. of the American Statistical Association , 93 , 1203 - 1215 ( 1998 )

17. Higuchi , T. : Self-organizing Time Series Model . In Sequential Monte Carlo Methods in Practice , Springer-Verlag New York, 429 - 444 ( 2001 )

18. Kitagawa , G. and Sato , S. : Monte Carlo Smoothing and Self-Organising Statespace Model . In Sequential Monte Carlo Methods in Practice , Springer-Verlag New York, 177 - 195 ( 2001 )

19. Kitagawa , G.: Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models . Journal of Computational and Graphical Statistics , 5 ( 1 ), 1 - 25 ( 1996 )

20. Heiner , M. , Lehrack , S. , Gilbert , D. , Marwan , W. : Extended Stochastic Petri Nets for Model-based Design of Wetlab Experiments . Edited by Priami, C. , et al., Trans. on Comput. Syst. Biol. XI, LNBI , 5750 , 138 - 163 ( 2009 )

21. Donaldson , R. , Gilbert , D. : A Monte Carlo Model Checker for Probabilistic LTL with Numerical Constraints , Technical report , University of Glasgow, Department of Computing Science ( 2008 )

22. Hu ¨rzeler and Hans , M. , Ku¨nsch, R.: Approximating and maximising the likelihood for a general state-space model . In: Sequential Monte Carlo Methods in Practice , Springer-Verlag New York, 159 - 175 ( 2001 )