Pacific Symposium on Biocomputing 14:227-238 (2009)

September 23, 2008

8:56

Proceedings Trim Size: 9in x 6in

psb09

PARAMETER ESTIMATION OF IN SILICO BIOLOGICALPATHWAYS WITH PARTICLE FILTERING TOWARDS APETASCALE COMPUTINGKAZUYUKI NAKAMURA1,∗, RYO YOSHIDA1, MASAO NAGASAKI2,
SATORU MIYANO2 AND TOMOYUKI HIGUCHI1
1

2

The Institute of Statistical Mathematics,
4-6-7 Minami-Azabu,
Minato-ku, Tokyo 1068569, JapanThe Institute of Medical Science, The University of Tokyo,
4-6-1 Shirokanedai,
Minato-ku, Tokyo 1088639, JapanThe aim of this paper is to demonstrate the potential power of large-scale particle filtering for the parameter estimations of in silico biological pathways wheretime course measurements of biochemical reactions are observable. The methodof 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-degreeof accuracy. In this paper, we address some implementation issues on large-scaleparticle filtering, and then, indicate the importance of large-scale computing forparameter learning of in silico biological pathways. We have tested the ability ofthe particle filtering with 108 Monte Carlo samples on the transcription circuitof 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. Introduction
Mathematical modeling of biological pathways has broad utility to understanding of complex networks of biochemical reactions in living cells. In
systems biology, in silico modeling of biological pathways has proven to be∗ to whom correspondence should be addressed: nakakazu@ism.ac.jp

Pacific Symposium on Biocomputing 14:227-238 (2009)

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 diﬀerential1,5,6 or diﬀerence2 – 4,7,8 equations that describe kinetic mechanism between biochemical reactants, e.g.
mRNA and protein. In order to conduct simulation experiments, it is essential to ﬁnd a set of model parameters and initial concentrations of reactants
such that the solution ﬁts observations. For example, one may ﬁnd eﬀective
values of rate constant parameters for which time course microarray data
are available.
More recently, several approaches have been proposed as a techniqueof estimating parameters of in silico pathway, e.g. genetic algorithm,8 unscented Kalman ﬁlter,6 and gradient methods.5 We also have been working
on this topic based on the method of particle ﬁltering,10 – 13 which has been
a popular technique for the problem of system identiﬁcations in the ﬁeldof statistical science and signal processing. While these existing methods
could cope with the estimation of relatively small number of parameters,
their ability is limited by the case that model contains much larger numberof unknown parameters.
As an illustration, consider the use of particle ﬁltering that approximates
the posterior distributions of model 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 togenerate an exceedingly large number of Monte Carlo samples in order tocover the overall parameter space. Such a problem, often called curse ofdimensionality, is also common in the existing methods.5,6,8 For example,
genetic algorithm explores values of parameters by alternating cross-overand 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 anexceedingly large amount of Monte Carlo samples, e.g. of order 108 or
more, for improving ability of particle ﬁltering. In our previous works,14,15
we applied the particle ﬁltering with 10, 000 Monte Carlo samples to theparameter estimation of transcription network of circadian clock. 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 ﬁltering has been extended
greatly so that 45 unknown parameters could be estimated with a high degree of accuracy. Particularly, the large-scale particle ﬁltering could reveal
clearly the shape of the posterior distributions over the 45 parameters. In

Pacific Symposium on Biocomputing 14:227-238 (2009)

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-modality and higher-ordermoments.
One more attractive feature of the particle ﬁltering is a high degreeof compatibility with parallel computing. This favorable feature increases
the potential of particle ﬁltering such that we can use a larger number ofparameters 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 ﬁrst give a brief explanation about the method of theparticle ﬁltering, and then, address some computational problems that arise
from the implementation of large-scale particle ﬁltering. The developmentof the 10 petaﬂops supercomputer has been going on as a national project.
It mainly targets high performance computing of life science simulators,
including molecular dynamics, biomechanics and biochemical simulations.
This work is a pilot study towards the estimation of large-scale in silico
models by the aid of petascale computing.
In Section 2, we show the method of particle ﬁltering after introducing
the state space model that relates in silico pathway models to time course
measurements. We also discuss some computational issues on the largescale particle ﬁltering. Section 3 shows the numerical experiments whichtest the applicability of the proposed approach on the analysis of transcription network of circadian clock. Finally, discussion and conclusion are given
in Section 4.
2. Method
2.1. Nonlinear State Space Model
Suppose that we model 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 vectorof model 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 methods are available to be associated with (1),for example, discretization in time. 9

Pacific Symposium on Biocomputing 14:227-238 (2009)

September 23, 2008

8:56

Proceedings Trim Size: 9in x 6in

psb09

To deal with an inevitable modeling error, we may rather use thestochastic diﬀerence equations, referred to as system model,
xt = f (xt−1 , θ) + v t ,(2)

where the noise component v t is added to (1). In what follows, we use thegeneral form (2) which contains (1) as a special case.
The measurement model, referred to as observation model, associates
the simulation model to time course data:y t = h(xt−1 ) + w t .

(3)Here, h(·) represents the function of measurement mechanism, which deﬁnes the conversion from unobserved concentration level xt to observed
time course data y t . wt denotes additive measurement noise. Initial statevariables x0 are assumed to follow a certain distribution p(x0 ). The probabilistic model consisting of (2) and (3) is referred to as the nonlinear state
space model.10 – 12
Of interest is then to ﬁnd values of parameters, θ, and unobserved
concentrations of reactants, xt , based on time course data y t , t = 1, · · · , T .
This paper focuses on the Bayesian estimations of model parameters, which
are based on the joint posterior distribution of xt , t = 0, · · · , T and θ conditional on the observed data. However, the Bayesian estimations usuallyrequire us to approximate the posterior distributions because their closed
forms are usually unavailable due to nonlinearity of models f .13,16
2.2. Particle Filtering
The particle ﬁltering is a sequential Monte Carlo technique 10 – 12 whichapproximates a joint posterior distribution of state vectors and parametersp(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 model parameters, respectively. Here, sets of state vectors and observations up to time t are denoted respectively by x0:t = {x0 , x1 , . . . , xt } andy 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 procedure of the particle ﬁltering is summarized as follows:
(i)

(i)

• Set t ← 1. Generate N Monte Carlo samples {x0|0 , θ0 }N
i=1 from
the prior distributions of initial state variables p(x0 ) and parameter
vector p(θ).

Pacific Symposium on Biocomputing 14:227-238 (2009)

September 23, 2008

8:56

Proceedings Trim Size: 9in x 6in

psb09• () Perform the following procedure:(i)– (Prediction step) After drawing {v t }N
i=1 , generate N predic(i)(i)
tion samples of state variables {xt|t−1 , θt−1 }N
i=1 according to
(i)

(i)

(i)

(i)

xt|t−1 = f (xt−1|t−1 , θt−1 ) + v t for i = 1, · · · , N .
– () (Filtering step) Execute the following procedure:(i)(i) Compute likelihood 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 ) forN (i)
i = 1, · · · , N such that i=1 qt = 1.
(i)(i)(iii) Generate {xt|t , θ t }Nby the resampling of
i=1
(i)

(i)

(i)

N
{xt|t−1 , θt−1 }N
i=1 with the probabilities {qt }i=1 .

• If t is equal to T , stop the algorithm. Otherwise, set t ← t + 1 andgo to ().
Here, the likelihood function p(y t |xt ) is deﬁned by the probability density
function, e.g. normal distribution, corresponding to the observation model(i)(3). The ﬁnal draws of {θT }N
i=1 generated from this procedure are approximately distributed according to the augmented posterior distributionp(x0:T , θ|y 1:T ). Indeed, it can be proved that the empirical distributions(i)given by N realization {θT }N
i=1 approach to the true distribution p(θ|y 1:T )with probability one as N goes to inﬁnity.13
As the number of unknown parameters becomes large, the method requires to generate an exceedingly large number of Monte Carlo samples inorder to cover the overall parameter space. For example, the transcription
regulation model of circadian clock 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
i=1 drawn from p(θ), the number of distinct
values in the resampled particles can diminish drastically as the time stepincreases. Hence, in order to estimate a large number of parameters, it isessential 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 timecomplexity of the ﬁltering step is O(N log N ) because the resampling algorithm requires O(log N ) evaluations of weight table to yield one ﬁlteredparticle (see Appendix A for details). We employ this type of resampling
procedure in the numerical experiment. Rather than performing the ran-

Pacific Symposium on Biocomputing 14:227-238 (2009)

September 23, 2008

8:56

Proceedings Trim Size: 9in x 6in

psb09

dom resampling algorithms, the ﬁltering step can be replaced by the resampling algorithm with the proportional-allocation that yields the ﬁlteredsamples by increasing and decreasing the number of prediction samples proportional to their weights. In that case, the time complexity can be reduced
to linear order.
Here, we brieﬂy discuss the low coupling property of the particle ﬁltering. In the prediction step, each predictive sample is updated independently. In the likelihood calculation (i) at the ﬁltering step, each likelihood
can be also computed independently. As a result, the tasks in each step
can be divided into parallel processors. On the other hand, the processes ofnormalization (ii) and resampling (iii) require additional parallel computing
technique in order to manage the communication among N Monte Carlosamples. In the step (ii), it is necessary that the likelihood values of N
samples are collected into a root 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 processors are available,
which indicates scalability of the particle ﬁltering in parallel computing.
3. Numerical Experiments
3.1. Implementation
Circadian clocks of living cells are controlled by oscillating feedback loop in
the transcription circuits of clock-control genes.2,19 The model that we consider was constructed by incorporating the regulatory mechanisms among
the ﬁve mRNAs (per, cry, rev/erb, clock, 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 consistsof the 12 ﬁrst-order stochastic diﬀerence 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 ,Pacific Symposium on Biocomputing 14:227-238 (2009)

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, clock, 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 constantparameters k· , the 7 threshold parameters s· , and the initial concentrationsof 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 clocks based on the time course gene expression data which weregenerated with Aﬀymetrix 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 ﬁve mRNAs (per, cry, rev/erb, clock, 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 allocated to the 700 simulation time steps by the conductionof time alignment. Note that in the current model (5), there are no regulatory mechanisms which yield oscillating concentrations of clock mRNA. As
a possible solution, we decided to treat the measurements of clock mRNA
concentrations as an input sequence. To this end, we interpolated the observed 12 data points of clock by using spline functions, and then, theinterpolated values were substituted into the corresponding state variables
which remained to be ﬁxed during the parameter estimation process.
Figure 1 shows the observed time course of gene expression proﬁles andthe result of simulation which was generated according to (5) with the predetermined parameters without inclusion of system noises vi,t s. This ﬁgureshows inconsistency in periods, phases and amplitude among the simulation
paths and observed proﬁles. The objective of this numerical experimentis to capture the observed oscillating patterns of circadian rhythm by theapplication of the particle ﬁltering where the currently obtained parameters
have the bias in the cyclic behavior.
Before proceeding to the particle ﬁltering, the observation noise wt was
set as the normal distribution N (0, 0.64) for per mRNA and N (0, 0.25)

Pacific Symposium on Biocomputing 14:227-238 (2009)

September 23, 2008

8:56

Proceedings Trim Size: 9in x 6in

psb09

PER/CRY

CLOCK

1.5
1.0

3.5

0.0

4.5

1.0

5.5

2.0

2.0

6.5

per

0

100

300

500

700

0

100

300

500

700

500

700

300

500

700

0

100

300

500

700

500

700

BMAL

1.0

0.5

1.5

1.0

2.0

1.5

2.5

2.0
1.5
1.0

700

1.0
100

REV/ERB

0.5

500

0.5
0

cry

2.0

300

300

1.5

1.5
0.5
100

100

bmal

1.0

2.5
1.5
0.5
0

0

rev/erb

PER

0

100

300

500

700

0

100

CRY

300

500

700

0

100

clock
1.2
0.0

0.5

0.4

1.0

0.8

1.5

2.0
1.5
1.0

300

BMAL/CLOCK

0

100

300

500

700

0

100

300

500

700

0

100

300

500

700

Figure 1. Result of initial simulation which was generated according to (5) with thepre-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 thesimulation shown in Figure 1 as the prior mean of the parameter θ andthe 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 , coeﬃcients k· and thresholdvalues s· were set to 0.01, 0.001 and 0.05, respectively.
In order to realize the implementation of the particle ﬁltering with8
10 samples, all particles were divided into 20 sets of 5 × 106 particles,
and allocated into 20 computational units. After repeating the particle ﬁltering with 5 × 106 samples for each computational unit, which
6(i,k)
yields {θT }10
i=1 for k = 1, · · · , 20, we can compute the estimator by
20
5×106 (i,k)θ̂ = 1018 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 theposterior distribution of the parameters kdi (i = 1, . . . , 12), where the sizesof 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 eﬃcient estimates by using such degenerated particles. On thePacific Symposium on Biocomputing 14:227-238 (2009)

Proceedings Trim Size: 9in x 6in

0.6

0.8

100000100000100000

50000

50000

0

0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

rev/erb

0.6

0.8

1.0

0.8

1.0

0.8

1.0

bmal

0.8

100000100000

50000

50000

0

0
0.0

1.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

cry

0.6

0.8

0.0

1.0

0.2

0.6

150000150000100000
50000

50000

100000
50000

0.4

BMAL

1000003.0e+072.0e+07
1.5e+071.0e+07

0.4

REV/ERB

0
0.0
0.2

0.4

0.6

0.8

0.4

0.6

0.8

1.0

0.0

0.2

CRY

3.0e+072.5e+072.0e+072.5e+07

1.5e+072.0e+071.0e+07

1.5e+07

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

1.0

BMAL/CLOCK

clock

BMAL/CLOCK
1500000.0

100000

1.0

50000

0.8

0

0

5.0e+06
0.6

0.0e+00
0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

0.2

0.4

0.6

0.8

0

5.0e+06

1.0e+071.0

0.6

0.0e+00
0.4

5.0e+06
0.8

1.0

1000000.4

2.5e+072.5e+072.0e+070.2

0.0e+00
0.6

0.8

50000
0.2

BMAL

3.0e+073.0e+072.5e+072.0e+07
1.5e+071.0e+075.0e+06

0.4

0.0

clock

0.0e+00

0.2

50000

1.0

1.5e+07
0.0

CRY

0.0

0.6

1500003.0e+072.0e+070.8

1.0e+071.0

0.4

1.5e+070.6

5.0e+06
0.8

0.2

PER

1.0e+070.4

0.0e+00

0.6

0.0

5.0e+06
0.2

3.0e+073.0e+072.5e+072.0e+07
1.5e+071.0e+075.0e+06

0.4

1.0

0.0e+00
0.0

REV/ERB

0.0e+00

0.2

0.8

1500001.0

0.6

100000

0.8

cry

0.0

0.4

50000

0.6

0.2

2.5e+072.5e+072.0e+07
1.5e+071.0e+075.0e+06
0.0e+00
0.4

0.0

bmal

3.0e+073.0e+072.5e+072.0e+07
1.5e+071.0e+075.0e+06

0.2

0

1.0

1500000.41500000.2

rev/erb

0.0e+00
0.0

1500003.0e+072.5e+072.0e+07
1.5e+071.0e+075.0e+06
0.0

1.0

CLOCK

0

0.8

0

0.6

PER

100000

0.4

50000

0.2

PER/CRY

0

0.0

per

0.0e+00

0.0e+00

0.0e+00

5.0e+06

5.0e+06

1.0e+071.0e+07

1.5e+07

1.5e+072.0e+072.0e+072.5e+072.5e+073.0e+073.0e+07

CLOCK

150000PER/CRY

150000per

psb09

1500008:56

150000

September 23, 2008

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

1.0

Figure 2. Approximated posterior distributions of the 12 “degradation rate” kd· , whichwere generated by using 108 (left) and 105 (right) particles.

other hand, the ability of the particle ﬁltering was improved greatly by theuse 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 averagesof 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 ﬁnal 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 periodicity of circadian rhythm whereas the diﬀerence in periodicity has been reduced greatly
by using 108 particles compared with the initially-chosen parameters (priormeans).
3.3. Estimation of computational time
We evaluate computational time for execution of large-scale particle ﬁltering
on upcoming peta-scale supercomputers.
We measured and estimated CPU times in the following three cases:

Pacific Symposium on Biocomputing 14:227-238 (2009)

Proceedings Trim Size: 9in x 6in

700

0

100

300

500

700

500

700

0

100

300

500

700

0

100

300

500

700

500

700

1.5
0.5

700

500

700

500

700

1.0
0.6
500

700

700

100

300

500

700

1.5
1.0

1.4
300

500

700

0.5
0

100

300

500

700

300

1.2
0

clock

1.0
100

100

0.8

2.5
500

0.6
0

0

BMAL

1.5
300

1.8

1.4
1.0
0.6
0.2
300

100

CRY

1.5

100

0

BMAL/CLOCK

1.0
0

300

3.5

1.5
1.0
0.5
700

0.5

Figure 3.

700

500

1.6

2.0
2.0
500

100

REV/ERB

1.5
300

0

cry

1.0
100

2.0

500

300

0.4

300

BMAL

1.5

300

100

0

100

300

BMAL/CLOCK

1.4

100

clock

1.0

100

0

bmal

1.0
0

2.5
0

CRY

0

700

1.0

700

1.5
700

500

2
500

1.0
500

300

1
300

2.0

1.5
1.0

300

100

rev/erb

0.8
100

REV/ERB

0.5

100

0

PER

0.4
0

cry

0

1.0

5.5
4.5
3.5

0.0
700

0.6

700

500

0.2

500

300

1.2

1.4
0.6
0.2
300

100

bmal

1.0

2.5
1.5

100

0

1.8

500

1.4

300

rev/erb

1.4

100

0.5
0

1.0

1.5
0

0.6

700

0.2

500

CLOCK

4.5

300

PER

1.6

100

PER/CRY

per

1.0

3.5

0.0
0

CLOCK

0.5

4.5

1.0

5.5

2.0

PER/CRY

2.0

per

psb09

3

8:56

2.0

September 23, 2008

0

100

300

500

700

0

100

300

500

700

State transition estimated by using 108 (left) and 105 (right) particles.

(a) 108 particles are divided into 20 sets of 5 million particles, andthen, repeat the particle ﬁltering 20 times on Opteron 2200 (about
5 gigaﬂops).
(b) Particle ﬁltering with 108 (100 million) samples is performed on one
petaﬂops machine.
(c) Particle ﬁltering with 109 (one billion) samples is performed on one
petaﬂops 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 petaﬂops 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.

Ave.
Min.
Max.

Opteron, 5mil. × 20

1 petaflops, 100 mil.

1 petaflops, 1 billion

6.8 × 105
3.7 × 105
2.0 × 1059050
2.7 × 102

1.0 × 103
5.6 × 102
3.0 × 103

Pacific Symposium on Biocomputing 14:227-238 (2009)

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
ﬁltering towards exhaustive explorations for kinetic parameters of in silico
biological pathways. The approach could oﬀer an overall view of posterior
distributions deﬁned over the 45 unknown parameters within a reasonable
computing time. To our knowledge, the use of one hundred million particles
is the ﬁrst attempt not only in systems biology but also in the other scientiﬁcﬁelds. Because the idea is too obvious, such a direct and straightforward
approach may have been overlooked 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 diﬃculty, we have to develop a methodology forincorporating the existing knowledge on biochemistry and systems biology.
According to the evaluation of computational time, it has been foundthat we can execute the particle ﬁltering even with one billion particles
within 20 minutes by the aid of one petaﬂops computer that will be accessible in Japan. The results of the experiments indicate that the large scaleparticle ﬁltering can be a promising approach to parameter estimation ofin silico biological pathways in the near future.
Appendix A. Resampling Algorithm
The ﬁltering step () of the algorithm that has O(N log N ) time complexity
is achieved as follows:(i)• Compute likelihood 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 stepis given by O(N ).
i
(i)
(j)
• Compute cumulative weight table q̂t for i = 0, · · · , N by j=1 q̂t(0)(q̂t = 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 ≤ q̂t by binary– Search the number k where q̂tsearch. Time complexity is O(N log N ).
(j)
(j)
(k)
(k)
– Set (xt|t , θt ) ← (xt|t−1 , θt−1 ), which requires O(N ).

Pacific Symposium on Biocomputing 14:227-238 (2009)

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 withmuch larger number of particles, we can replace the binary search algorithm
by the resampling with the proportional-allocation that yields the ﬁlteredsamples 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 Cell., 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’Alché-Buc, Bioinformatics, 23, 3209-3216(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. Thiagarajan, Bioinformatics, 22, e271-e280 (2006).
9. B.-E. Perrin, L. Ralaivola, A. Mazurie, S. Bottani, J. Mallet, F. d’Alché-Buc,
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, Parallel Programming with MPI, Morgan Kaufmann (1996).
18. Y. Fujii, Y. Okitsu, H. Matsuno, S. Miyano and S. T. Inoue, A new regulatoryinteractions 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).