September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 PARAMETER ESTIMATION OF IN SILICO BIOLOGICAL PATHWAYS WITH PARTICLE FILTERING TOWARDS A PETASCALE COMPUTING KAZUYUKI NAKAMURA1,, RYO YOSHIDA1, MASAO NAGASAKI2, SATORU MIYANO2 AND TOMOYUKI HIGUCHI1 1 The Institute of Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 1068569, Japan 2 The Institute of Medical Science, The University of Tokyo, 4-6-1 Shirokanedai, Minato-ku, Tokyo 1088639, Japan The aim of this paper is to demonstrate the potential power of large-scale particle filtering for the parameter estimations of in silico biological pathways where time course measurements of biochemical reactions are observable. The method of particle filtering has been a popular technique in the field of statistical science, which approximates posterior distributions of model parameters of dynamic system by using sequentially-generated Monte Carlo samples. In order to apply the particle filtering to system identifications of biological pathways, it is often needed to explore the posterior distributions which are defined over an exceedingly high-dimensional parameter space. It is then essential to use a fairly large amount of Monte Carlo samples to obtain an approximation with a high-degree of accuracy. In this paper, we address some implementation issues on large-scale particle filtering, and then, indicate the importance of large-scale computing for parameter learning of in silico biological pathways. We have tested the ability of the particle filtering with 108 Monte Carlo samples on the transcription circuit of circadian clock that contains 45 unknown kinetic parameters. The proposed approach could reveal clearly the shape of the posterior distributions over the 45 dimensional parameter space. 1. Intro duction Mathematical mo deling of biological pathways has broad utility to understanding of complex networks of bio chemical reactions in living cells. In systems biology, in silico mo deling of biological pathways has proven to be to whom correspondence should be addressed: nakakazu@ism.ac.jp September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 a promising approach with successful applications.1 ­ 4 Usually, a pathway model is comprised of a set of differential1,5,6 or difference2 ­ 4,7,8 equations that describe kinetic mechanism between bio chemical reactants, e.g. mRNA and protein. In order to conduct simulation experiments, it is essential to find a set of mo del parameters and initial concentrations of reactants such that the solution fits observations. For example, one may find effective values of rate constant parameters for which time course microarray data are available. More recently, several approaches have been proposed as a technique of estimating parameters of in silico pathway, e.g. genetic algorithm,8 unscented Kalman filter,6 and gradient metho ds.5 We also have been working on this topic based on the metho d of particle filtering,10 ­ 13 which has been a popular technique for the problem of system identifications in the field of statistical science and signal pro cessing. While these existing metho ds could cope with the estimation of relatively small number of parameters, their ability is limited by the case that mo del contains much larger number of unknown parameters. As an illustration, consider the use of particle filtering that approximates the posterior distributions of mo del parameters based on a set of Monte Carlo samples. The size of random numbers usually lies in 103 to 105 .13 As the number of unknown parameters becomes large, it is required to generate an exceedingly large number of Monte Carlo samples in order to cover the overall parameter space. Such a problem, often called curse of dimensionality, is also common in the existing metho ds.5,6,8 For example, genetic algorithm explores values of parameters by alternating cross-over and mutation of initially-generated seeds for unknown parameters.8 In its applications, it is a key to success whether or not an initial set of candidate parameters lies in the region close to the optimal values. To address such a limitation, we adopt the direct way that uses an exceedingly large amount of Monte Carlo samples, e.g. of order 108 or more, for improving ability of particle filtering. In our previous works,14,15 we applied the particle filtering with 10, 000 Monte Carlo samples to the parameter estimation of transcription network of circadian clo ck. Unfortunately, the number of parameters that were estimated was limited to only a few. However, as will be demonstrated, by increasing the number of samples to 108 particles, the ability of the particle filtering has been extended greatly so that 45 unknown parameters could be estimated with a high degree of accuracy. Particularly, the large-scale particle filtering could reveal clearly the shape of the posterior distributions over the 45 parameters. In September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 doing so, we aim to obtain an overall view of high-dimensional parameter space of in silico biological pathways, e.g. multi-mo dality and higher-order moments. One more attractive feature of the particle filtering is a high degree of compatibility with parallel computing. This favorable feature increases the potential of particle filtering such that we can use a larger number of parameters to deal with higher dimensional parameter space, e.g. of order 102 , by exploiting the parallel computation schemes and computational resources which will be available in the near future. In this paper, we first give a brief explanation about the metho d of the particle filtering, and then, address some computational problems that arise from the implementation of large-scale particle filtering. The development of the 10 petaflops supercomputer has been going on as a national pro ject. It mainly targets high performance computing of life science simulators, including molecular dynamics, biomechanics and bio chemical simulations. This work is a pilot study towards the estimation of large-scale in silico mo dels by the aid of petascale computing. In Section 2, we show the metho d of particle filtering after intro ducing the state space mo del that relates in silico pathway mo dels to time course measurements. We also discuss some computational issues on the largescale particle filtering. Section 3 shows the numerical experiments which test the applicability of the proposed approach on the analysis of transcription network of circadian clo ck. Finally, discussion and conclusion are given in Section 4. 2. Metho d 2.1. Nonlinear State Space Model Suppose that we mo del successive change in concentration levels of biochemical reactants by a set of rate equations xt = f (xt-1 , ) , (1) where the p dimensional vector xt contains concentration levels of all variables, e.g. mRNAs and proteins, at time step t and denotes a vector of mo del parameters. The state variable xt which is assumed to be unobservable evolves during successive discrete time points from t - 1 to t according to the p dimensional vector-valued function f . For ordinary differential equations, several metho ds are available to be asso ciated with (1), for example, discretization in time. 9 September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 To deal with an inevitable mo deling error, we may rather use the sto chastic difference equations, referred to as system mo del, xt = f (xt-1 , ) + v t , (2) where the noise component v t is added to (1). In what follows, we use the general form (2) which contains (1) as a special case. The measurement mo del, referred to as observation mo del, asso ciates the simulation mo del to time course data: y t = h(xt-1 ) + w t . (3) Here, h(·) represents the function of measurement mechanism, which defines the conversion from unobserved concentration level xt to observed time course data y t . wt denotes additive measurement noise. Initial state variables x0 are assumed to follow a certain distribution p(x0 ). The probabilistic mo del consisting of (2) and (3) is referred to as the nonlinear state space mo del.10 ­ 12 Of interest is then to find values of parameters, , and unobserved concentrations of reactants, xt , based on time course data y t , t = 1, · · · , T . This paper fo cuses on the Bayesian estimations of mo del parameters, which are based on the joint posterior distribution of xt , t = 0, · · · , T and conditional on the observed data. However, the Bayesian estimations usually require us to approximate the posterior distributions because their closed forms are usually unavailable due to nonlinearity of mo dels f .13,16 2.2. Particle Filtering The particle filtering is a sequential Monte Carlo technique 10 ­ 12 which approximates a joint posterior distribution of state vectors and parameters p(x0:T , |y 1:T ) p(y 1:T |x1:T )p(x1:T | , x0 )p(x0 )p( ) , (4) where p(x0 ) and p() denote the prior distributions of initial concentrations and mo del parameters, respectively. Here, sets of state vectors and observations up to time t are denoted respectively by x0:t = {x0 , x1 , . . . , xt } and y 1:t = {y1 , y 2 , . . . , y t }. Denote N Monte Carlo samples from the conditional distribution (i) (i) p(xt , |y 1:k ) by (xt|k , k ), i = 1, · · · , N . Then, the pro cedure of the particle filtering is summarized as follows: · Set t 1. Generate N Monte Carlo samples {x0|0 , 0 }N 1 from i= the prior distributions of initial state variables p(x0 ) and parameter vector p( ). (i) (i) September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 · ( ) Perform the following pro cedure: ­ (Prediction step) After drawing {v t }N 1 , generate N predici= (i) (i) tion samples of state variables {xt|t-1 , t-1 }N 1 according to i= xt|t-1 = f (xt-1|t-1 , t-1 ) + v t for i = 1, · · · , N . ­ ( ) (Filtering step) Execute the following pro cedure: (i) Compute likeliho o d function of each particle p(y t |xt|t-1 ), i = 1, · · · , N . (i) (i) (ii) Calculate the normalizing weights qt p(y t |xt|t-1 ) for N (i) i = 1, · · · , N such that i=1 qt = 1. (i) (i) (iii) Generate {xt|t , t }N 1 by the resampling of i= {xt|t-1 , t-1 }N 1 with the probabilities {qt }N 1 . i= i= · If t is equal to T , stop the algorithm. Otherwise, set t t + 1 and go to ( ). Here, the likelihood function p(y t |xt ) is defined by the probability density function, e.g. normal distribution, corresponding to the observation mo del (i) (3). The final draws of {T }N 1 generated from this pro cedure are api= proximately distributed according to the augmented posterior distribution p(x0:T , |y 1:T ). Indeed, it can be proved that the empirical distributions (i) given by N realization {T }N 1 approach to the true distribution p(|y 1:T ) i= with probability one as N go es to infinity.13 As the number of unknown parameters becomes large, the metho d requires to generate an exceedingly large number of Monte Carlo samples in order to cover the overall parameter space. For example, the transcription regulation mo del of circadian clo ck that will be shown in Section 3 contains 45 unknown parameters in total. Particularly, the problem of "particle degeneracy" is caused by a small number of Monte Carlo samples (Chapter 14 of Ref. 13). This is due to the fact that during the repeated resampling (i) of initially generated {0 }N 1 drawn from p(), the number of distinct i= values in the resampled particles can diminish drastically as the time step increases. Hence, in order to estimate a large number of parameters, it is essential to use a fairly large amount of Monte Carlo samples. Both the time complexity and space complexity of the prediction step are linear order of the number of particles. On the other hand, the time complexity of the filtering step is O(N log N ) because the resampling algorithm requires O(log N ) evaluations of weight table to yield one filtered particle (see Appendix A for details). We employ this type of resampling pro cedure in the numerical experiment. Rather than performing the ran(i) (i) (i) (i) (i) (i) (i) (i) (i) September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 dom resampling algorithms, the filtering step can be replaced by the resampling algorithm with the proportional-allo cation that yields the filtered samples by increasing and decreasing the number of prediction samples proportional to their weights. In that case, the time complexity can be reduced to linear o rder. Here, we briefly discuss the low coupling property of the particle filtering. In the prediction step, each predictive sample is updated independently. In the likeliho o d calculation (i) at the filtering step, each likeliho o d can be also computed independently. As a result, the tasks in each step can be divided into parallel pro cessors. On the other hand, the pro cesses of normalization (ii) and resampling (iii) require additional parallel computing technique in order to manage the communication among N Monte Carlo samples. In the step (ii), it is necessary that the likeliho o d values of N samples are collected into a ro ot CPU. To implement a parallel computing for such interacted tasks, we can use the divide and conquer algorithm.17 By the aid of it, time complexity of normalization and resampling can be reduced up to O((N/M ) log M ) where M CPU pro cessors are available, which indicates scalability of the particle filtering in parallel computing. 3. Numerical Exp eriments 3.1. Implementation Circadian clo cks of living cells are controlled by oscillating feedback lo op in the transcription circuits of clo ck-control genes.2,19 The mo del that we consider was constructed by incorporating the regulatory mechanisms among the five mRNAs (per, cry, rev/erb, clo ck, bmal), the translated proteins (PER, CRY, REV/ERB, CLOCK, BMAL) and the two protein complexes (PER/CRY, CLOCK/BMAL), which were reported by Ref. 18. It consists of the 12 first-order sto chastic difference equations as, x1,t+1 = x1,t - kd1 x1,t + ktc1 I (x5,t < s1,5 )I (x12,t > s1,12 ) + kc1 + v1,t+1 , x2,t+1 = x2,t - kd2 x2,t + ktl2 x1,t - kb2,4 x2,t x4,t + v2,t+1 , x3,t+1 = x3,t - kd3 x3,t + ktc3 I (x5,t < s3,5 )I (x12,t > s3,12 ) + kc3 + v3,t+1 , x4,t+1 = x4,t - kd4 x4,t + ktl4 x3,t - kb2,4 x2,t x4,t + v4,t+1 , x5,t+1 = x5,t - kd5 x5,t + kb2,4 x2,t x4,t + v5,t+1 , x6,t+1 = x6,t - kd6 x6,t + ktc6 I (x5,t < s6,5 ) + kc6 + v6,t+1 , x7,t+1 = x7,t - kd7 x7,t + ktl7 x6,t + v7,t+1 , x8,t+1 = x8,t - kd8 x8,t + kc8 + v8,t+1 , September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 x9,t+1 = x9,t - kd9 x9,t + ktl9 x8,t - kb9,11 x9,t x11,t + v9,t+1 , x10,t+1 = x10,t - kd10 x10,t+1 + ktc10 I (x5,t > s10,5 )I (x7,t < s10,7 ) +kc10 + v10,t+1 , x11,t+1 = x11,t - kd11 x11,t + ktl11 x10,t - kb9,11 x9,t x11,t + v11,t+1 , x12,t+1 = x12,t - kd12 x12,t + kb9,11 x9,t x11,t + v12,t+1 , (5) in order corresponding to per, PER, cry, CRY, PER/CRY, rev/erb, REV/ERB, clo ck, CLOCK, bmal, BMAL, and CLOCK/BMAL, respectively. Here, vi,t , i = 1, · · · , 12, stand for the Gaussian noises with mean 0 and variance 2.0 × 10-7. Of interest is then to estimate the 26 rate constant parameters k· , the 7 threshold parameters s· , and the initial concentrations of 12 reactants xi,0 , i = 1, · · · , 12. Note that we did not estimate values of kd8 and kc8 for the reason shown below. Ueda et al.19 studied the transcription regulations of mammalian circadian clo cks based on the time course gene expression data which were generated with Affymetrix mouse high-density oligonucleotide probe array. Time course data were measured at equally-spaced 12 time points during 44 hours (about two days). We extracted gene expression measurements of the five mRNAs (per, cry, rev/erb, clo ck, bmal) from Figure 1a and the supplementary data (http://sirius.cdb.riken.jp/ MouseLiver/MouseLiverCCG(020707).html 19 ), where the set of 12 time points were allo cated to the 700 simulation time steps by the conduction of time alignment. Note that in the current mo del (5), there are no regulatory mechanisms which yield oscillating concentrations of clo ck mRNA. As a possible solution, we decided to treat the measurements of clo ck mRNA concentrations as an input sequence. To this end, we interpolated the observed 12 data points of clo ck by using spline functions, and then, the interpolated values were substituted into the corresponding state variables which remained to be fixed during the parameter estimation pro cess. Figure 1 shows the observed time course of gene expression profiles and the result of simulation which was generated according to (5) with the predetermined parameters without inclusion of system noises vi,t s. This figure shows inconsistency in perio ds, phases and amplitude among the simulation paths and observed profiles. The ob jective of this numerical experiment is to capture the observed oscillating patterns of circadian rhythm by the application of the particle filtering where the currently obtained parameters have the bias in the cyclic behavior. Before pro ceeding to the particle filtering, the observation noise wt was set as the normal distribution N (0, 0.64) for per mRNA and N (0, 0.25) September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 per 6.5 2.0 PER/CRY 2.0 CLOCK 5.5 1.0 4.5 0 .0 0 100 300 500 700 3.5 0 100 300 500 700 1.0 0 1.5 100 300 500 700 PER 1.5 2.5 rev/erb 1.5 bmal 1.5 1.0 0.5 0 2.0 100 300 500 700 0.5 0 100 300 500 700 2.0 0.5 0 1.0 100 300 500 700 cry 2.5 REV/ERB BMAL 1 .5 1 .0 2.0 0.5 1.5 0 100 300 500 700 1.0 0 100 300 500 700 0.5 0 1.0 1.5 100 300 500 700 CRY 2.0 1.5 clock 1.2 BMAL/CLOCK 1 .5 1.0 1.0 0.5 0 100 300 500 700 0 100 300 500 700 0.0 0 0.4 0.8 100 300 500 700 Figure 1. Result of initial simulation which was generated according to (5) with the pre-determined parameters without inclusion of system noises vi,t s. The observed 12 data points are indicated by dots. for the other reactants. We adopted the values that are employed in the simulation shown in Figure 1 as the prior mean of the parameter and the initial state x0 . Because the parameters and the states should be nonnegative, the prior distributions were set to truncated normal distributions. The prior variances of initial concentration x·,0 , co efficients k· and threshold values s· were set to 0.01, 0.001 and 0.05, respectively. In order to realize the implementation of the particle filtering with 8 10 samples, all particles were divided into 20 sets of 5 × 106 particles, and allo cated into 20 computational units. After repeating the particle filtering with 5 × 106 samples for each computational unit, which (i,k) 06 yields {T }1=1 for k = 1, · · · , 20, we can compute the estimator by i 20 5×106 (i,k) 1 ^ = 108 k=1 i=1 T . Due to the independence of the computational ^ units, the integrated estimator still preserves the statistical consistency. 3.2. Results Figure 2 shows the histograms of Monte Carlo samples drawn from the posterior distribution of the parameters kdi (i = 1, . . . , 12), where the sizes of Monte Carlo samples were set to 108 and 105 , respectively. It is observed from the right panel of Figure 2 that the number of distinct values declines drastically for the estimation with 105 particles. Obviously, we would not obtain any efficient estimates by using such degenerated particles. On the September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 per 3.0e+07 3.0e+07 2.5e+07 2.5e+07 PER/CRY 3.0e+07 2.5e+07 CLOCK 150000 per 150000 PER/CRY 150000 CLOCK 2.0e+07 2.0e+07 2.0e+07 100000 1.5e+07 1.5e+07 1.0e+07 1.0e+07 1.0e+07 1.5e+07 100000 50000 5.0e+06 5.0e+06 0.0e+00 0.0e+00 0.0e+00 5.0e+06 0 50000 0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0 0.0 50000 100000 0.2 0.4 0.6 0.8 1.0 PER 3.0e+07 3.0e+07 2.5e+07 2.5e+07 rev/erb 3.0e+07 bmal 2.5e+07 150000 PER 150000 rev/erb 150000 bmal 2.0e+07 2.0e+07 1.5e+07 1.5e+07 1.0e+07 1.0e+07 1.5e+07 2.0e+07 100000 100000 1.0e+07 5.0e+06 5.0e+06 0.0e+00 0.0e+00 0.0e+00 5.0e+06 50000 50000 0 0 0 0.0 50000 100000 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 cry 3.0e+07 3.0e+07 2.5e+07 REV/ERB 3.0e+07 150000 150000 2.5e+07 2.0e+07 2.5e+07 100000 100000 2.0e+07 1.5e+07 1.5e+07 1.0e+07 1.5e+07 2.0e+07 50000 1.0e+07 50000 5.0e+06 5.0e+06 0.0e+00 0.0e+00 5.0e+06 1.0e+07 0 0.0e+00 0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0 50000 100000 150000 BMAL cry REV/ERB BMAL 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 CRY 3.0e+07 3.0e+07 clock 3.0e+07 150000 2.5e+07 2.5e+07 2.5e+07 150000 100000 2.0e+07 2.0e+07 1.5e+07 1.5e+07 1.5e+07 2.0e+07 100000 50000 1.0e+07 1.0e+07 5.0e+06 5.0e+06 5.0e+06 1.0e+07 50000 0.0e+00 0.0e+00 0.0e+00 0 0 0 50000 100000 150000 BMAL/CLOCK CRY clock BMAL/CLOCK 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 2. Approximated posterior distributions of the 12 "degradation rate" kd· , which were generated by using 108 (left) and 105 (right) particles. other hand, the ability of the particle filtering was improved greatly by the use of 108 particles as the histograms of the Monte Carlo samples reveal the shapes of posterior distributions clearly as shown in the left panel of Figure 2. Figure 3 shows the results of simulations based on the estimated parameters and initial state variables, which were computed as the averages of 108 and 105 particles, respectively. In the case of 108 particles, simulation paths were generated from the estimated parameters in each of the 20 computational units, and then, the generated 20 paths were averaged in order to compute the final estimates shown in Figure 3. The observed data points are also shown in the Figure 3. The simulated concentrations corresponding to 105 particles preceded the observed perio dicity of circadian rhythm whereas the difference in perio dicity has been reduced greatly by using 108 particles compared with the initially-chosen parameters (prior means). 3.3. Estimation of computational time We evaluate computational time for execution of large-scale particle filtering on upcoming peta-scale supercomputers. We measured and estimated CPU times in the following three cases: September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 2.0 2.0 per PER/CRY CLOCK 2 .0 per 5.5 PER/CRY CLOCK 1.5 5.5 1.5 1 .0 1.0 1.0 4.5 4.5 0.5 0.0 0.0 0 100 300 500 700 3.5 0 100 300 500 700 1.6 0 100 300 500 700 0 100 300 500 700 1.4 3.5 0 100 300 500 700 1.8 0.5 0 1.0 100 300 500 700 PER 2 .5 1.4 rev/erb bmal PER 3 rev/erb bmal 1.2 1 .5 1.0 0.8 2 1.0 0.6 0.5 0.4 1 0.6 0 100 300 500 700 2.5 0.2 0 100 300 500 700 0 100 300 500 700 2 .0 0 100 300 500 700 4.5 0.2 0 100 300 500 700 0.6 0 1.0 1.4 100 300 500 700 cry 1 .5 REV/ERB 2.0 BMAL cry REV/ERB 1.6 BMAL 2.0 1 .5 3.5 1.0 1.5 2.5 1.5 0 .5 1.0 1.0 1.5 0 2.0 100 300 500 700 1.0 0 100 300 500 700 1.4 0 100 300 500 700 0 .5 0 100 300 500 700 0 100 300 500 700 1.4 0.4 0 0.8 1.2 100 300 500 700 CRY 1.5 clock BMAL/CLOCK 1.8 CRY 1.5 clock BMAL/CLOCK 1 .5 1.0 1.0 1.4 1.0 0.6 1 .0 0.5 1 .0 0.5 0.2 0 100 300 500 700 0 100 300 500 700 0 100 300 500 700 0.6 0 100 300 500 700 0 100 300 500 700 0.2 0 0.6 1.0 100 300 500 700 Figure 3. State transition estimated by using 108 (left) and 105 (right) particles. (a) 108 particles are divided into 20 sets of 5 million particles, and then, repeat the particle filtering 20 times on Opteron 2200 (about 5 gigaflops). (b) Particle filtering with 108 (100 million) samples is performed on one petaflops machine. (c) Particle filtering with 109 (one billion) samples is performed on one petaflops machine. For the evaluations of (b) and (c), the CPU times were estimated based on the observed CPU time of the experiment (a). The results of the tests are summarized in Table 1. The current implementation that divides 108 particles into 20 computational subunits requires about eight days on average, which can be reduced to 90 seconds by the aid of one petaflops computer. Even with one billion particles, the estimated CPU time can be within 20 minutes. Table 1. Observed and estimated CPU times (second). The left column shows the measured CPU time on Opteron 2200 and 20 sets of 5 million particles. The middle and right columns show the estimated CPU times for 100 million and one billion particles. Opteron, 5mil. × 20 Ave. Min. Max. 6.8 × 105 3.7 × 105 2.0 × 105 1 petaflops, 100 mil. 90 50 2.7 × 102 1 petaflops, 1 billion 1.0 × 103 5.6 × 102 3.0 × 103 September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 4. Conclusion This paper has shown the promising approach based on large scale particle filtering towards exhaustive explorations for kinetic parameters of in silico biological pathways. The approach could offer an overall view of posterior distributions defined over the 45 unknown parameters within a reasonable computing time. To our knowledge, the use of one hundred million particles is the first attempt not only in systems biology but also in the other scientific fields. Because the idea is to o obvious, such a direct and straightforward approach may have been overlo oked so far. A limitation of our approach is the sensitivity of the estimation results to the choice of the prior distributions, that is, if the prior distributions are fairly inappropriate, we would not obtain reasonable posterior estimations. In order to overcome such a difficulty, we have to develop a metho dology for incorporating the existing knowledge on bio chemistry and systems biology. According to the evaluation of computational time, it has been found that we can execute the particle filtering even with one billion particles within 20 minutes by the aid of one petaflops computer that will be accessible in Japan. The results of the experiments indicate that the large scale particle filtering can be a promising approach to parameter estimation of in silico biological pathways in the near future. App endix A. Resampling Algorithm The filtering step ( ) of the algorithm that has O(N log N ) time complexity is achieved as follows: · Compute likeliho od function of each particle p(y t |xt|t-1 ), i = 1, · · · , N . Time complexity of this step is given by O(N ). (i) (i) · Calculate the normalizing weights qt = p(y t |xt|t-1 ) for i = N (i) 1, · · · , N such that i=1 qt = 1. Time complexity of this step is given by O(N ). i (i) (j ) ^ · Compute cumulative weight table qt for i = 0, · · · , N by j =1 qt ^ (qt = 0). Time complexity is given by O(N ). ^ · For each j = 1, · · · , N ­ Generate random number u from uniform distribution U (0, 1]. Time complexity is O(N ). (k-1) (k ) < u qt by binary ^ ­ Search the number k where qt ^ search. Time complexity is O(N log N ). (j ) (j ) (k ) (k ) ­ Set (xt|t , t ) (xt|t-1 , t-1 ), which requires O(N ). (0) (i) September 23, 2008 8:56 Proceedings Trim Size: 9in x 6in psb09 Only the binary searching requires O(N log N ). In order to deal with much larger number of particles, we can replace the binary search algorithm by the resampling with the proportional-allo cation that yields the filtered samples by increasing and decreasing the number of prediction samples proportional to their weights. References 1. K. C. Chen, L. Calzone, A. Csikasz-Nagy, F. R. Cross, B. Novak and J. J. Tyson, Mol Biol Cel l., 15, 3841-3862 (2004). 2. H. Matsuno, R. Murakami, R. Yamane, N. Yamasaki, S. Fujita, H. Yoshimori and S. Miyano, Pac Symp Biocomput., 8, 152-163 (2003). 3. A. Doi, S. Fujita, H. Matsuno, M. Nagasaki and S. Miyano, In Silico Biol., 4, 271-291 (2004). 4. A. Doi, M. Nagasaki, H. Matsuno and S. Miyano, In Silico Biol., 6, 1-13 (2006). 5. K. Fujarewicz, M. Kimmel, T. Lipniacki and A. Swierniak, IEEE/ACM Transaction on Computational Biology and Bioinformatics, 4, 322-335 (2007). 6. M. Quach, N. Brunel and F. d'Alch´-Buc, Bioinformatics, 23, 3209-3216 e (2007). 7. I. Nachman, A. Regev and N. Friedman, Bioinformatics, 20, i248-i256 (2004). 8. G. Koh, H. F. C. Teong, M-V ment, D. Hsu and P.S. Thiagara jan, Bioinformatics, 22, e271-e280 (2006). 9. B.-E. Perrin, L. Ralaivola, A. Mazurie, S. Bottani, J. Mallet, F. d'Alch´-Buc, e Bioinformatics, 19, ii386-ii349 (2003). 10. G. Kitagawa, J. Computational and Graphical Stat. 5, 1-25 (1996). 11. N. J. Gordon, D. J. Salmond and A. F. M. Smith IEE Proceedings-F, 140, 107-113 (1993). 12. T. Higuchi and G. Kitagawa, IEICE Tras. Inf. Syst., E83-D, 36-43 (2000). 13. A. Doucet and N. De Freitas, Springer-Verlag (2001). 14. M. Nagasaki, R. Yamaguchi, R. Yoshida, S. Imoto, A.Doi, Y.Tamada, H. Matsuno, S. Miyano and T. Higuchi, Genome Inform. 17(1), 46-61 (2006). 15. S. Tasaki, M. Nagasaki, M. Oyama, H. Hata, K. Ueno, R. Yoshida, T. Higuchi, S. Sugano, S. Miyano, Genome Inform. 17(2), 226-238 (2006). 16. G. Kitagawa, J. Amer. Stat. Assoc., 93, 1203-1215 (1998). 17. Peter Pacheco, Paral lel Programming with MPI, Morgan Kaufmann (1996). 18. Y. Fujii, Y. Okitsu, H. Matsuno, S. Miyano and S. T. Inoue, A new regulatory interactions suggested by simulations for circadian genetic control mechanism in mammals, Genome Inform., available at http://www.jsbi.org/journal/ GIW0, (2004). 19. H. R. Ueda, W. Chen, A. Adachi, H. Wakamatsu, S. Hayashi, T. Takasugi, M. Nagano, K. Nakahama, Y. Suzuki, S. Sugano, M. Iino, Y. Shigeyoshi and S. Hashimoto, Nature, 418, 534-539 (2002).