Bayesian Parameter Inference and Model Selection by Population Annealing in Systems Biology

15 Pages·2014·0.91 MB·English

Bayesian Parameter Inference and Model Selection by
Population Annealing in Systems Biology
Yohei Murakami*
DepartmentofBiophysics,DivisionofBiology,GraduateSchoolofScience,KyotoUniversity,Kyoto,Japan
Abstract
Parameter inference and model selection are very important for mathematical modeling in systems biology. Bayesian
statistics can be used to conduct both parameter inference and model selection. Especially, the framework named
approximateBayesiancomputationisoftenusedforparameterinferenceandmodelselectioninsystemsbiology.However,
MonteCarlomethodsneedstobeusedtocomputeBayesianposteriordistributions.Inaddition,theposteriordistributions
ofparametersaresometimesalmostuniformorverysimilartotheirpriordistributions.Insuchcases,itisdifficulttochoose
one specific value of parameter with high credibility as the representative value of the distribution. To overcome the
problems, we introduced one of the population Monte Carlo algorithms, population annealing. Although population
annealingisusuallyusedinstatisticalmechanics,weshowedthatpopulationannealingcanbeusedtocomputeBayesian
posterior distributions in the approximate Bayesian computation framework. To deal with un-identifiability of the
representativevaluesofparameters,weproposedtorunthesimulationswiththeparameterensemblesampledfromthe
posterior distribution, named ‘‘posterior parameter ensemble’’. We showed that population annealing is an efficient and
convenientalgorithmtogenerateposteriorparameterensemble.Wealsoshowedthatthesimulationswiththeposterior
parameterensemblecan,notonlyreproducethedatausedforparameterinference,butalsocaptureandpredictthedata
which was not used for parameter inference. Lastly, we introduced the marginal likelihood in the approximate Bayesian
computation framework for Bayesian model selection. We showed that population annealing enables us to compute the
marginallikelihoodintheapproximateBayesiancomputationframeworkandconductmodelselectiondependingonthe
Bayes factor.
Citation:MurakamiY(2014)BayesianParameterInferenceandModelSelectionbyPopulationAnnealinginSystemsBiology.PLoSONE9(8):e104057.doi:10.
1371/journal.pone.0104057
Editor:MariaAnisimova,SwissFederalInstituteofTechnology(ETHZurich),Switzerland
ReceivedMarch24,2014;AcceptedJuly8,2014;PublishedAugust4,2014
Copyright:(cid:2)2014YoheiMurakami.Thisisanopen-accessarticledistributedunderthetermsoftheCreativeCommonsAttributionLicense,whichpermits
unrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalauthorandsourcearecredited.
DataAvailability:Theauthorconfirmsthatalldataunderlyingthefindingsarefullyavailablewithoutrestriction.Partofrelevantdataarewithinthepaperand
itsSupportingInformationfiles.OtherrelevantdataareavailablefromthewebpageofthegroupofNoguchietal.(2013)Mol.Syst.Biol.study.
Funding:ThisworkwaspartlysupportedbyaGrant-in-AidforScientificResearchonInnovativeAreasfromtheMinistryofEducation,Culture,Sports,Science,
andTechnology(MEXT)ofJapan,byGrantsforExcellentGraduateSchools,MEXT,JapanandbyGrant-in-AidforScientificResearchoftheMEXT.Thefundershad
noroleinstudydesign,datacollectionandanalysis,decisiontopublish,orpreparationofthemanuscript.
CompetingInterests:Theauthorhasdeclaredthatnocompetinginterestsexist.
*Email:[email protected]
Introduction AIC [9] has already used to compare a number of mathematical
models insystems biology[10–12].
Mathematicalmodelingisaveryusefulandpowerfulapproach
Alternativetothosemethodsandapproaches,Bayesianstatistics
in systems biology [1,2]. Mathematical models used in systems
enables us to conduct both model selection [11,13–19] and
biology are often represented by ordinary or partial differential
parameter inference [13,15,16,20,21] under the coherent philos-
equations. These differential equations contain a number of
ophy. For Bayesian model selection and parameter inference in
parameters which represent the rates of biochemical reactions or
systems biology, the framework named approximate Bayesian
amounts of components (proteins, mRNAs etc). An appropriate computation (ABC) or likelihood-free computation is often used
mathematical model together with appropriate values of param- [14,15,18,22–27]. ABC is very useful when likelihood functions
eters enables us to explain, understand and predict the target are analytically or computationally intractable, or it takes much
biological phenomena in system level. However, we may have a cost for evaluation [14,15,23,24]. Instead of evaluating the
number of competing and potential mathematical models to likelihood function, ABC compares the observed data and the
explain the observed experimental data. In addition, concrete simulated data, andgives us theapproximated posterior distribu-
values of parameters in mathematical models are often not well- tion [14,15].
knowninpreviousexperimentalliteratures.Insuchcases,weneed However, to compute Bayesian posterior distributions, Monte
toconductmodelselectionandparameterinferencebysomesort Carlo methods often need to be used. As Bayesian parameter
of systematic procedures. inferenceandmodelselectionalgorithmsintheABCframework,
Forparameterinference,wecanuselikelihoodbasedapproach Markov chain Monte Carlo (MCMC) [15,22,24] and sequential
as a statistical approach [3]. In addition, many optimization MonteCarlo(SMC)[15,18,23]havealreadyused.Itisknownthat
algorithms have already been developed and used to conduct the efficiency of ABC-MCMC algorithm reduces when the
parameterinferenceinsystemsbiology[4–8].Formodelselection, sampler is trapped in a low probability area [23]. To overcome
PLOSONE | www.plosone.org 1 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
the problem, ABC-SMC algorithm was developed [15,18,23]. As the weighting function, so-called indicator function is often
However, even though the algorithms work successfully, the used[14,22–24,26]. Indicator functionisrepresented as
posterior distributions of parameters may be almost uniform or
very similar to their prior distributions [13]. In such cases, it is f ðD jD ,h,MÞ!I(cid:2)d(cid:2)D ,D (cid:3)ƒe,M(cid:3):
difficult to choose one specific value of parameter with high w obs sim obs simjh
credibilityastherepresentativevalueoftheposteriordistribution.
d(D ,D )isadistancemeasurebetweentheobserveddataand
obs sim|h
To overcome the problems, we introduced one of the Monte the simulated data. e§0 is a tolerance. The indicator function
Carlo algorithms, population annealing [28,29]. Population
equals to 1 if the observed data and the simulated data are close
annealing is a so-called population Monte Carlo algorithm [30] (!e) and 0 if not (.e). If e is small enough, the augmented
assameasSMC.Althoughpopulationannealingisusuallyusedin
posteriordistributionisagoodapproximationofthetrueposterior
statisticalmechanics[28,31],weshowedthatpopulationannealing
distribution [14,24]. Thus, we can use the augmented posterior
can be used to compute Bayesian posterior distributions for
distribution alternative to the true posterior distribution for
parameter inference and model selection in the ABC framework.
parameter inference.
To deal with un-identifiability of the representative values of
parameters,insteadofchoosingonespecificvalueofparameterin
Bayesian model selection
the posterior distribution, we ran the simulations with the
The marginal likelihood P(D |M) plays an important role in
parameter ensemble sampled from the posterior distribution, obs
Bayesianmodelselection[11,17,37,38].Themarginallikelihoodis
named ‘‘posterior parameter ensemble’’. We propose this
represented as
approach is valid if our purpose of parameter inference is to
reproduce or predict the system dynamics, not to estimate the
ð
correct values of parameters. We showed that population PðD jMÞ~ fðD jh,MÞpðhjMÞdh:
annealing is an efficient and convenient algorithm to generate obs obs
H
theposteriorparameterensemble.Inaddition,weshowedthatthe
simulations with the posterior parameter ensemble can, not only InBayesianmodelselection,amodelMisavariable.Tocompare
reproducethedatausedforparameterinference,butalsocapture modelsM1andM2,wecanusetheBayesfactor[37]represented
andpredictthedatawhichwasnotusedforparameterinference. as
Lastly, we introduced the marginal likelihood in the ABC
framework for Bayesian model selection. We showed that PðD jM Þ
population annealing enables us to compute the marginal B12~PðDobsjM1Þ:
obs 2
likelihoodintheABCframeworkandtoconductmodelselection
depending on the Bayes factor. The validity of our propositions FromtheBayes’theorem,theBayesfactoristheratioofposterior
wasfirstlyjudgedbyapplyingourmethodtothefeed-forwardloop oddsand prior oddsasfollows[37]:
networkmotifmodels[32–36],secondlytotheinsulindependent
AKT pathway model [47,48].
PðD jM Þ PðM jD Þ=PðM jD Þ
B ~ obs 1 ~ 1 obs 2 obs :
12 PðD jM Þ PðM Þ=PðM Þ
Methods obs 2 1 2
P(M) is the prior probability of the model M. P(M|D ) is the
Bayesian parameter inference obs
posteriorprobabilityofthemodelM.Whenthepriorprobabilities
For Bayesian parameter inference, under the given likelihood
of the competing models are equal (P(M )=P(M )), the Bayes
function f(D |h,M) and the prior distribution of parameters 1 2
obs factor equalstotheratio oftheposterior probabilities
p(h|M), we try to obtain the posterior distribution of parameters
p(h|D ,M),represented as
obs
PðD jM Þ PðM jD Þ
B ~ obs 1 ~ 1 obs :
12 PðD jM Þ PðM jD Þ
pðhjD ,MÞ obs 2 2 obs
obs
fðD jh,MÞpðhjMÞ Depending on the Bayes factor, we can conduct Bayesian model
~ obs !fðD jh,MÞpðhjMÞ:
PðD jMÞ obs selection [15,17,37,38].
obs
To conduct Bayesian model selection in the ABC framework,
M is a model and fixed in parameter inference. h is a set of weneedtodefinethemarginallikelihoodintheABCframework.
parameters (i.e. multidimensional vector) in the model M. D is For definition, westart fromthe augmented posterior probability
obs
an observed experimental data. In the equation, P(D |M) is a of a model represented as
obs
normalization constant, and alsocalled asa marginal likelihood.
IntheABCframework[14,24,26],insteadofthetrueposterior
P ðM,D jD Þ
ABC sim obs
distribution,wetrytoobtaintheaugmentedposteriordistribution
Ð
pABC(h,Dsim|Dobs,M) [24]. The augmented posterior distribution ~PðMÞ HfwðDobsjDsim,h,MÞfðDsimjh,MÞpðhjMÞdh:
can berepresented as PðD Þ
obs
On another front, wehavetherelationship
p ðh,D jD ,MÞ!f ðD jD ,h,MÞfðD jh,MÞpðhjMÞ:
ABC sim obs w obs sim sim
D is the simulated data, also called as the auxiliary parameter PðMjD Þ PðD jMÞ
sim obs ~ obs :
[24]. fw(Dobs|Dsim,h,M) is the weighting function [24]. The PðMÞ PðDobsÞ
weighting function is set to be large value when the observed
dataandthesimulateddataare‘‘close’’,smallvaluewhentheyare Thenumeratoroftherighthandsideoftheequationcorresponds
‘‘distant’’,andconstantwhentheyare‘‘equal’’(D =D )[24]. totheconventionalmarginallikelihood.Wecanobtainthesimilar
obs sim
PLOSONE | www.plosone.org 2 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
relationship with theaugmented posterior distribution as annealing,wedefinetheintermediatedistributionsf (x)assame
IM
as annealed importance sampling [39]. The first intermediate
PABCðM,DsimjDobsÞ sdaismtrpibliuntgioinsefaIsMy0.(Tx)hieslsaesttitnotetrhmeepdrioabteabdiilsittryibduitsitorinbuftioNn(xo)fiswsehtictho
PðMÞ IM
the target distribution i.e. Bayesian posterior distribution in
Ð
f ðD jD ,h,MÞfðD jh,MÞpðhjMÞdh Bayesianapproach.Inpopulationannealing,thefirstintermediate
~ H w obs sim sim :
PðD Þ distribution is gradually changed to the last intermediate
obs
distributionthroughanumberofintermediatedistributionsf n(x)
IM
Thus, we can assume the marginal likelihood in the ABC (0,n,N)betweenthefirstone(n=0)andthelastone(n=N).In
framework as thecaseofstatisticalmechanics,theintermediatedistributionsare
canonical distributions with different temperatures of the system.
ð Thus,annealingisdirectlycorrespondstothegradualdecreaseof
fwðDobsjDsim,h,MÞfðDsimjh,MÞpðhjMÞdh: temperature of canonical distributions. Application to Bayesian
H statistics mimics the process [39]. In population annealing,
Thevalidityoftheapproximationcanbeconfirmedbysettingthe particles are moved and their weights are changed to follow the
weighting functiontotheindicatorfunction as intermediate distribution in each annealing step. Population
annealing algorithm proceeds asfollows[28,29]:
PA1.Generatexk,f 0(x)(1!k!K)independentlyandsetthe
fwðDobsjDsim,h,MÞ~IðDobs~DsimÞ: initialvalueofweighttoIMw k=1/K,whereKisthetotalnumberof
0
particles.
The indicator function equals to 1 if the observed data equals to
PA2. Repeat the following procedure from n=1 to n=N for
the simulated data and 0 if not. By this setting, the marginal
each particle independently.
likelihood in the ABC framework is consistent with the conven-
tional marginal likelihoodas follows: (a) Update weightsby following equation
ð
fwðDobsjDsim,h,MÞfðDsimjh,MÞpðhjMÞdh logwk~logwk zlog fInM(cid:2)xk(cid:3) (cid:4)wk~wk fInM(cid:2)xk(cid:3)(cid:5),
H n n{1 fn{1ðxkÞ n n{1fn{1ðxkÞ
ð IM IM
~ IðDobs~DsimÞfðDsimjh,MÞpðhjMÞdh andnormalize.
H
ð (b) Update xk independently by finite number of MCMC
~ fðDobsjh,MÞpðhjMÞdh movements as the stationary distribution is consistent with
H then-thintermediatedistribution f n(x).
IM
~PðDobsjMÞ (c) Resample particles at appropriate timing by following
procedure.
Thus,ifwesettheweightingfunctiontotheindicatorfunctionas
For eachparticle, set theprobability as
fwðDobsjDsim,h,MÞ!I(cid:2)d(cid:2)Dobs,Dsimjh(cid:3)ƒe,M(cid:3), wk
Pk~ n :
and e is small enough, the marginal likelihood in the ABC XK wk
k~1 n
frameworkisagoodapproximationoftheconventionalmarginal
likelihood. This is same as the case of the augmented posterior Sample new particles with replacement following the
distribution for parameter inference in the ABC framework. probabilities andset theweight aswnk=1/K.
Therefore, the Bayes factor can be defined and calculated in the (d) Set n=n+1 andreturn to(a).
ABC framework as
InthePAalgorithm,wetransposedtheMCMCprocessandthe
Ð resamplingprocessintheoriginalPAalgorithm[28].AsMCMC
f ðD jD ,h,M ÞfðD jh,M ÞpðhjM Þdh
BA12BC~ÐHHfwwðDoobbssjDssiimm,h,M12ÞfðDssiimmjh,M12ÞpðhjM12Þdh: p(TroexcetssS1in) [P4A0,241(]b.),TwoeuusseeptohpeuMlateiotrnopaonlins-eHalainstgingins athlgeorAithBmC
framework,weneedtouseABC-MCMC(SeeTextS1)inPA2(b)
DependingontheBayesfactorBABC ,wecanconductBayesian
12 [22,24]. To decide the appropriate timing of resampling in PA2
model selection in theABCframework.
(c), we used effective sample size (ESS) as same as Sisson et al’s
study [23].ESS isdefined as
Population annealing: algorithm
Population annealing [28,29] is one of the population Monte
1
Carlo algorithm [30] as same as SMC. Thus, in principle, ESS~ :
populationannealingcanavoiditssamplerbeingtrappedinalow XK (cid:2)wk(cid:3)2
k~1
probability area as same as SMC [23]. Although population
annealingisusuallyusedinstatisticalmechanicstosamplefroma ESS represents the bias of particles’ weights. Particles are
canonicaldistribution[28,31],weshowthatpopulationannealing resampled when the ESS value falls below a threshold. In this
can beusedtocomputeBayesian posterior distributions. study,wesetthethresholdtoK/2.TheoutputofthePAalgorithm
Population annealing uses particles to approximately represent istheKparticleswhicheachofthemhasxkassociatedwithweight
a target probability distribution. Each particle has one set of w k (1!k!K). These weights are set to follow the target
N
concrete values of a multidimensional variable x, which corre- distribution f N(x). By counting xk associated with weights w k,
IM N
sponds to a sample from the target distribution. For population wecanobtain thetarget distribution f N(x).
IM
PLOSONE | www.plosone.org 3 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
Populationannealing:applicationtoBayesianparameter of resampling, the number of non-zero weight particles before
inference resamplingkresampleshouldbememorized.Bymultiplyingkresample/
K tothe number of particles at theend of the algorithm, wecan
Population annealing can be used to conduct Bayesian
estimate the marginal likelihood in the ABC framework. In this
parameter inference. In the case of parameter inference, a
study, we adopted this solution. This solution is based on the
multidimensional variable x corresponds to a set of parameters h
similar idea of the ‘‘ABC-SMC approximation of the marginal
and a fixed model M. The target distribution for sampling is the
likelihood’’ [18]. Another solution against resampling is to
posteriordistributionofparametersp(h|D ,M)(SeeTextS2asa
obs
eliminate the resampling process from population annealing. In
examplebasedon[17,39])ortheaugmentedposteriordistribution
this study, we did not adopt that solution. This is because the
p (h,D |D ,M) in the ABC framework. In the ABC
ABC sim obs
resampling process is considered to be important to avoid the
framework, we firstly set the weighting function to the indicator
fluctuation of the weight factor and conduct stable calculations
function. Then we can define the intermediate distribution as
[28].
follows:
The computation of the number of particles satisfying the
condition d(D ,D )!e can be also done by ABC rejection
fInMðhÞ!I(cid:2)d(cid:2)Dobs,Dsimjh(cid:3)ƒen,M(cid:3)fðDsimjh,MÞpðhjMÞ sampler (TextobsS1)sim[15,23,24,42]. Fixing the total number of
ðe w:::we w:::we Þ samplingtrials,thenumberofacceptedparticlesinABCrejection
0 n N
sampler is proportional to the marginal likelihood in the ABC
As the tolerance e decreases, the intermediate distribution framework. However, it is known that the acceptance rate often
n
gradually changes from the first intermediate distribution (e=e ) gets lower in the case that the prior distribution is very different
0
to the last intermediate distribution (e=e ) corresponding to the fromtheposteriordistribution[15].Wecanavoidtheproblemby
N
augmented posterior distribution of parameters as the similar population annealing, because the intermediate distributions
mannerinABC-SMC [15,18].Theconcrete scheduleofe isset gradually changes from the prior distribution to the posterior
n
depending on theproblem. distribution inpopulation annealing.
Population annealing: application to Bayesian model Simulations with the posterior parameter ensemble
selection Inthisstudy,insteadofsamplinganumberofparametersfrom
the posterior distribution, we ran the simulations with all the
Population annealing can be used to conduct Bayesian model
output particles of population annealing. This is fundamentally
selection in the ABC framework. We firstly set the weighting
same as the simulations with the posterior parameter samples
functiontotheindicatorfunction.Underthissetting,withalarge
[11,13]. In population annealing, each particle at the end of the
numberK,themarginallikelihoodintheABCframeworkcanbe
algorithm (output particles of population annealing) has concrete
approximated as follows:
values of parameters hk associated with weight w k (1!k!K).
N
These weights are set to follow the target posterior distribution.
ð
fwðDobsjDsim,h,MÞfðDsimjh,MÞpðhjMÞdh Thus, the simulations with hk weighted by wNk for all the K
H particles correspond to the simulations with the parameter
ð ensemblesampledfromtheposteriordistribution.Practically,this
~ I(cid:2)d(cid:2)Dobs,Dsimjh(cid:3)ƒe,M(cid:3)fðDsimjh,MÞpðhjMÞdh can be done with use of the same particles used to obtain the
H
posterior distribution of parameters. In this manner, population
~EfðDsimjh,MÞpðhjMÞ(cid:6)I(cid:2)d(cid:2)Dobs,Dsimjh(cid:3)ƒe,M(cid:3)(cid:7) annealing is a convenient algorithm to generate the posterior
parameter ensemble andconduct theensemble simulation.
&1XK I(cid:2)d(cid:2)D ,D (cid:3)ƒe,M(cid:3)
K obs simjh,k Results
k~1
E[?] indicates an expectation value. D is a sample from Flow of the tests
sim|h,k
f(D |h,M)p(h|M).Thus,underafixedtotalnumberofparticles Inthissection,weexplaintheflowofthetestswhichstartfrom
sim
K,thenumberofparticlessatisfyingtheconditiond(D ,D )!e thenextsection.Totestthevalidityofourmethod,weappliedour
obs sim
isproportionaltothemarginallikelihoodintheABCframework. method to the two test problems. In the first test, we introduced
Populationannealingenablesustocountthenumberofparticles twomathematicalmodelsofthefeed-forwardloopnetworkmotif
satisfying the condition. Because we use the indicator function, [32–36]inthenext‘‘Mathematicalmodelsandsimulationsetting
eachparticlehasacommonvalueofweight(non-zeroweight)or of the feed-forward loop network motif’’ section. Then, we
zero weight in population annealing. Thus, starting with a fixed artificially generated the observed data from one of the models
total number of particles K, the number of particles which have with the given answer values of parameters. Next, in the
non-zeroweightattheendofpopulationannealingisproportional ‘‘Parameterinference,reproductionandpredictionofthedynam-
tothemarginallikelihoodintheABCframework.Thisisbecause ics of the feed-forward loop models’’ section, we firstly tried to
the particles at the end of the algorithm follow the last estimate the answer values of parameters by computing the
intermediate distribution corresponding tothetarget distribution posterior distribution of parameters. Secondly, we tested whether
thesimulationswithposteriorparameterensemblecanreproduce
the data used for parameter inference. In addition, we tested
fN ðhÞ!p ðh,D jD ,MÞ!IðdðD ,D Þƒe ,MÞ
IM ABC sim obs obs sim N whether the simulations with posterior parameter ensemble can
fðD jh,MÞpðhjMÞ: correctlypredictthenewlygenerateddatawhichwasnotusedfor
sim
parameter inference. Next, in the ‘‘Model selection between the
However, even though starting with a fixed total number of feed-forward loop models’’ section, we tested whether the true
particles, resampling processes recover the number of non-zero model that the observed data was generated from was correctly
weightparticlestothetotalnumberofparticlesK.Thus,eachtime selected or not by computing the Bayes factor. Lastly, in the
PLOSONE | www.plosone.org 4 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
remainingsections,weshowedtheresultsofthesecondtest.Inthe infer these seven values. In model selection, these seven
second test, we used the insulin dependent AKT pathway model parameters were free parameters. Remaining Hill coefficient was
[47,48].Inthetest,wetriedtoreproducetheexperimentaltime- set to h=2. This value was always given and fixed in this test
series data of phosphorylated AKT. We also conducted Bayesian study. With these answer values and the fixed value of Hill
model selection between the wild type AKT pathway model and coefficient, the coherent FFL and the incoherent FFL models
themutantAKTpathwaymodel.Importantly,weusedtheopen showeddistinctdynamicsofZinresponsetothestepstimulation
experimental datain thesecondtest. of X (Figure 1.A, B). Secondly, we added Gaussian noise with
mean 0 and variance 0.01 as observation noise to the time-series
Mathematicalmodelsandsimulationsettingofthefeed- data of Z in Figure 1.A. The generated data was represented as
forward loop network motif Dobs={Ztobs,t=1,10}(redpointsinFigure 1.C,concretevalues
Mathematical models of the feed-forward loop network are shown in Table S1). The observed data Dobs was used for all
motif. Forthefirsttest,weusedthemathematicalmodelsofthe the computations of parameter inference and model selection in
feed-forward loop network motif [32–36]. Network motifs are thefirst test withFFLmodels.
building blocks of transcription network found in diverse Prior distribution and weighting function. For Bayesian
organisms [34]. One of the significant network motifs is the parameter inference and model selection, we set the prior
feed-forwardloop(FFL),whichconsistsofthreecomponents,X,Y distribution of each parameter to independently follow the
and Z. There are totally eight possible structures of FFL, which uniformdistributiononacommonlogarithmicscale.Logarithmic
can be divided into coherent FFL or incoherent FFL [33,34]. In scale is often used in systems biology [13,14,43,44]. The upper-
addition,thereareANDlogicandORlogicfortheactivationofZ boundandthelower-boundofeachuniformdistributionwereset
byXandY[33,34].Inthecurrentteststudy,weusedoneofthe to 10-folds value and 1/10-folds value of each answer value of
coherent FFL and one of the incoherent FFL with AND logic parameter respectively i.e. log10aY , U[0.1,10], log10aZ ,
(Figure 1.A, B). Their dynamics can be represented by the U[0.1,10], log10bY , U[0.1,10], log10bZ , U[0.1,10], log10KXY
following ordinary differential equations [33]. ,U[0.01,1],log10KXZ,U[0.01,1],log10KYZ,U[0.05,5].Here,
Coherent FFL U represents theuniform distribution.
Inthisteststudy,wesettheweightingfunctiontotheindicator
function[14,22–24,26]represented as
(cid:4) X (cid:5)h
dY K
dt ~bY (cid:4)XYX (cid:5)h{aYY fwðDobsjDsim,h,MÞ!I(cid:2)d(cid:2)Dobs,Dsimjh(cid:3)ƒe,M(cid:3),
1z
KXY andset thedistance d(Dobs,Dsim|h)as
(cid:4) X (cid:5)h (cid:4) Y (cid:5)h d(cid:2)Dobs,Dsimjh(cid:3)~X10 (cid:8)Zotbs{Zstimjh(cid:9)2:
dZ K K t~1
~b XZ YZ {a Z
dt Z (cid:4) X (cid:5)h (cid:4) Y (cid:5)h Z D ={Zt , t=1,10} is a simulated time-series data of Z
1z 1z sim|h sim|h
K K with h.
XZ YZ
Numerical simulation. The total number of particles in
population annealing was set to K=100000. As the default
Incoherent FFL annealing schedule, the tolerance was gradually decreased as
follows:
(cid:4) X (cid:5)h
dY~b KXY {a Y ðe0,e1,e2,e3,e4,e5Þ~ð?, 1,0:5,0:25,0:2, 0:15Þ
dt Y (cid:4) X (cid:5)h Y
1z The proposal distribution of ABC-MCMC in population anneal-
K
XY ingwasset totheuniformdistribution ona commonlogarithmic
scale.Inconcreteterms,ateachstepofABC-MCMC,oneofthe
parameters was randomly chosen, and the uniform random
(cid:4) X (cid:5)h number between 20.25 to 0.25 was added on a common
dZ K 1 logarithmic scale. For each particle, ABC-MCMC movements in
~b XZ {a Z
dt Z (cid:4) X (cid:5)h (cid:4) Y (cid:5)h Z populationannealingweresetto7steps.Thisisthenumberofthe
1z 1z inferred parameters.
K K
XZ YZ For time-series calculations, the ordinary differential equations
werenumericallysolvedbythefourth-orderRunge-Kuttamethod
In the equations, X is an input signal, set to a step function or a
with a time step of 0.01. Both of the initial amounts of Y and Z
pulseinthetest.YandZarevariables.Eachmodelcontainseight
weresetto0inallthecalculations.Xwassettothestepfunction
parameters,a ,a ,b ,b ,K ,K ,K andh.Inthisteststudy,
Y Z Y Z XY XZ YZ (X=1 during the entire simulation) or the pulse (X=1 until
we did not specifically define the units of these parameters and
time=5,then X=0).
time forsimplicity.
Generationoftheartificialobserveddata. To use in the
test,weartificiallygeneratedtheobserveddatafromthecoherent Parameter inference, reproduction and prediction of the
FFL model (Figure 1.A). To generate the data, we firstly set the dynamics of the feed-forward loop models
answer values of parameters h =(a , a , b , b , K , K , Coherent FFL model. Firstly, we conducted parameter
answer Y Z Y Z XY XZ
K )=(1,1,1,1,0.1,0.1,0.5).Inparameterinference,wetriedto inference of the coherent FFL model by population annealing.
YZ
PLOSONE | www.plosone.org 5 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
Figure1.Feed-forwardloopnetworkmotifmodels.(A)StructureanddynamicsofthecoherentFFLmodelinresponsetothestepstimulation
ofX.(B)StructureanddynamicsoftheincoherentFFLmodelinresponsetothestepstimulationofX.(C)Generatedobserveddata.Theobserved
data(redpoints)wasgeneratedfromthecoherentFFLmodelbyaddingobservationnoiseintothegreen-coloredtrajectory.Thegreen-colored
trajectoryissameasthatofin(A).
doi:10.1371/journal.pone.0104057.g001
The computed joint posterior distribution of the inferred trajectories. Reproductions of the observed data used for
parameters was marginalized and shown in Figure2. In each parameter inference were shown in Figure 3. We ran the
marginal distribution, the red-colored class corresponds to the simulations in response to the step stimulation of X (Figure 3.A)
answervalueofeachparameter.InFigure 2,a andb seemedto as same as the generation process of the observed data. In
Z Z
be inferred well to some extent compared to other parameters. Figure 3.B, the area of the probability density consists of the
However,marginaldistributionsofotherparameterswerealmost simulated trajectories (blue-colored area) could capture the
similar to the uniform distributions, which are same as the prior observed data used for parameter inference (red points). When
distributionsinthisstudy.Ifourpurposeofparameterinferenceis theannealingschedulewaschangedtothesmallertolerances(e ,
0
the estimation of the correct values of unknown parameters with e , e , e , e , e )=(‘, 1, 0.5, 0.25, 0.15, 0.11), the area of the
1 2 3 4 5
high credibility, we failed in accomplishing our purpose in this probabilitydensitygotnarrower(Figure 3.C).Whentheannealing
case. schedule was changed to the larger tolerances (e , e , e , e , e ,
0 1 2 3 4
However,ifourpurposeisthereproductionorthepredictionof e )=(‘, 2, 1, 0.75, 0.5, 0.25), the area of the probability density
5
thesystemdynamics,insteadofchoosingonesetofrepresentative got broader and captured all the red points (Figure3.D). As
valuesofparameters,wecanrunthesimulationswiththeposterior shown,theareaoftheprobabilitydensitydiffersdependingonthe
parameterensemble.Inthisstudy,weranthesimulationswithuse annealing schedules. Thisresult isunsurprising because thelast e
of all the 100000 output particles (parameter sets) computed by values strongly restrict the acceptable trajectories to the observed
population annealing. The simulations with posterior parameter data. Thus, we are recommended to try a number of annealing
ensemble give the probability density consists of the simulated schedulestochecktheinfluenceoftheschedulesonthesimulated
PLOSONE | www.plosone.org 6 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
Figure2.MarginalprobabilitydistributionsoftheparametersinthecoherentFFLmodel.Thejointprobabilitydistributionapproximated
by100000particleswasmarginalizedforeachparameter.Red-coloredclassesrepresentthe‘‘answer’’valueofeachparameter.
doi:10.1371/journal.pone.0104057.g002
data. In either case, the posterior parameter ensemble could used for parameter inference. We newly generated the observed
reproduce theobserved datawell. data of Z (red points in Figure 4, concrete values are shown in
Furthermore, we examined whether the simulations with the Table S2) in response to the pulse stimulation of X (Figure 4.A)
posteriorparameterensemblecancapturethedatawhichwasnot adding the same Gaussian noise in Figure 1.C. Then we ran the
Figure 3. Reproduction of the observed data by the coherent FFL model. Simulations with the posterior parameter ensemble of the
coherentFFLmodelinresponsetothestepstimulationofX.(A)InputstepstimulationofX.(B)Annealingschedule:(e ,e,e,e ,e ,e)=(‘,1,0.5,
0 1 2 3 4 5
0.25,0.2,0.15).(C)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,1,0.5,0.25,0.15,0.11).(D)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,2,1,0.75,0.5,
0 1 2 3 4 5 0 1 2 3 4 5
0.25).Blue-coloredareaistheprobabilitydensityconsistsofthesimulatedtrajectories.Redpointsaretheobserveddata.
doi:10.1371/journal.pone.0104057.g003
PLOSONE | www.plosone.org 7 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
simulationsinresponsetothepulsestimulationofXwithuseofall values.ThisisoppositetothatofthecoherentFFL(Figure 2).b
Y
the 100000 particles which were the output of population and K also showed the opposite tendency between the
YZ
annealing. In Figure 4.B, the probability density consists of the incoherent FFL model and the coherent FFL model (Figure 2).
simulatedtrajectoriescouldcaptureandpredicttheobserveddata Thesedifferenceswereconsideredtocomefromthedifferenceof
whichwasnotusedforparameterinference.Inaddition,whenthe the interaction from Y to Z (Figure1.A, B). Although there were
annealing schedule was changed to the smaller tolerances (e , e , these kinds of small differences, most of the marginal posterior
0 1
e , e , e , e )=(‘, 1, 0.5, 0.25, 0.15, 0.11) (Figure4.C) or the distributions of the incoherent FFL model were similar to the
2 3 4 5
larger tolerances (e , e , e , e , e , e )=(‘, 2, 1, 0.75, 0.5, 0.25) uniform distributions. This result is similar to the result of the
0 1 2 3 4 5
(Figure 4.D), the areas of the probability densities could capture coherent FFLmodel.
and predicttheobserved datawith different rigors. In addition, as shown in Figure 6.B, C and D, the simulations
As a whole, the simulations with the posterior parameter with the posterior parameter ensemble of the incoherent FFL
ensemble could, not only reproduce the data used for parameter modelinresponsetothestepstimulationofX(Figure6.A)showed
inference, but also capture and predict the data which was not similarresultstothoseofthecoherentFFLmodel(Figure 3).We
used for parameter inference. Because the posterior parameter again note that the observed data used for parameter inference
ensemble corresponds to the output particles of population was generated from the coherent FFL model, not from the
annealing, population annealing is an efficient and convenient incoherent FFL model. This result indicates that, by setting the
algorithm for thesimulationswith posterior parameter ensemble. values of parameters suitably, even the false model (incoherent
Incoherent FFL model. Next, we conducted parameter FFL model) can reproduce the observed data with comparable
inference of the incoherent FFL model by population annealing. level to the true model (coherent FFL model). This result also
Note that parameter inference was conducted with the observed emphasizes the importance of the concrete values of kinetic
data generated from the coherent FFL model in Figure 1.C, not parameters for the dynamics of the system [45], not only the
fromtheincoherentFFLmodel.Thejointposteriordistributionof network structures. One attention point is the case of the
theinferredparameterswasmarginalizedandshowninFigure 5. annealing schedule (e , e , e , e , e , e )=(‘, 1, 0.5, 0.25, 0.15,
0 1 2 3 4 5
We can see the distinct tendencies of marginal distributions 0.11) (Figure6.C). In this case, the probability density became
betweenthecoherentFFLmodelandtheincoherentFFLmodel. slightlyparabolic.Thismightbeakindofover-fittingtothedata.
Forexample,themarginaldistributionofa hadatailinsmaller
Z
Figure4.PredictionoftheobserveddatabythecoherentFFLmodel.Simulationswiththeposteriorparameterensembleofthecoherent
FFLmodelinresponsetothepulsestimulationofX.(A)InputpulsestimulationofX.(B)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,1,0.5,0.25,0.2,
0 1 2 3 4 5
0.15).(C)Annealingschedule:(e ,e ,e ,e ,e ,e)=(‘,1,0.5,0.25,0.15,0.11).(D)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,2,1,0.75,0.5,0.25).Blue-
0 1 2 3 4 5 0 1 2 3 4 5
coloredareaistheprobabilitydensityconsistsofthesimulatedtrajectories.Redpointsaretheobserveddata.
doi:10.1371/journal.pone.0104057.g004
PLOSONE | www.plosone.org 8 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
Figure 5. Marginal probability distributions of the parameters in the incoherent FFL model. The joint probability distribution
approximatedby100000particleswasmarginalizedforeachparameter.
doi:10.1371/journal.pone.0104057.g005
However,thepredictionofZdynamicsinresponsetothepulse observeddatainresponsetothepulsestimulationofX,wemight
stimulation of X (Figure7.A) did not succeed well (Figure 7.B, C beabletonaivelyselectthetruemodel(coherentFFLmodel)ina
andD).AfterthedecreaseofXattime=5,partofthetrajectories visual way. However, it is not always possible to obtain the
ofZcouldnotcapturetheobserveddatawell.Thus,ifwehavethe
Figure6. Reproductionof theobserveddatabytheincoherentFFLmodel.Simulationswiththeposteriorparameterensembleofthe
incoherentFFLmodelinresponsetothestepstimulationofX.(A)InputstepstimulationofX.(B)Annealingschedule:(e ,e ,e ,e ,e ,e)=(‘,1,0.5,
0 1 2 3 4 5
0.25,0.2,0.15).(C)Annealingschedule:(e ,e ,e ,e ,e ,e)=(‘,1,0.5,0.25,0.15,0.11).(D)Annealingschedule:(e ,e ,e ,e ,e ,e)=(‘,2,1,0.75,0.5,
0 1 2 3 4 5 0 1 2 3 4 5
0.25).Blue-coloredareaistheprobabilitydensityconsistsofthesimulatedtrajectories.Redpointsaretheobserveddata.
doi:10.1371/journal.pone.0104057.g006
PLOSONE | www.plosone.org 9 August2014 | Volume 9 | Issue 8 | e104057
PopulationAnnealingforParameterInferenceandModelSelection
Figure 7. Prediction of the observed data by the incoherent FFL model. Simulations with the posterior parameter ensemble of the
incoherentFFLmodelinresponsetothepulsestimulationofX.(A)InputpulsestimulationofX.(B)Annealingschedule:(e ,e ,e ,e ,e ,e)=(‘,1,0.5,
0 1 2 3 4 5
0.25,0.2,0.15).(C)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,1,0.5,0.25,0.15,0.11).(D)Annealingschedule:(e ,e,e ,e ,e ,e)=(‘,2,1,0.75,0.5,
0 1 2 3 4 5 0 1 2 3 4 5
0.25).Blue-coloredareaistheprobabilitydensityconsistsofthesimulatedtrajectories.Redpointsaretheobserveddata.
doi:10.1371/journal.pone.0104057.g007
convenient data for model selection. To deal with this kind of value of the last e got smaller, BABC got larger. This result
CI
problem, wecanconduct Bayesian model selection. indicatesthat,asthemodelsneedtoreproducetheobserveddata
morerigorously,thecoherentFFLmodelisselectedwithstronger
Model selection between the feed-forward loop models evidence against the incoherent FFL model. In addition, BABC
CI
Next, we conducted Bayesian model selection, comparing the wasalwayslargerthan1.Thisresultisconsistentwiththefactthat
coherentFFLmodelwiththeincoherentFFLmodel.Toconduct theobserveddatausedformodelselectionwasgeneratedfromthe
Bayesian model selection, we computed the Bayes factor in the coherent FFL model, and the prediction that the coherent FFL
ABC framework represented as model must be selected in this case. Thus, Bayesian model
selection by populationannealing isvalid.
WealsoconductedBayesianmodelselectionbyABCrejection
BABC
CI sampler(TextS1).Forcomparison,thevalueofeinABCrejection
~ÐHÐIH(cid:2)Id(cid:2)(cid:2)dD(cid:2)oDbso,bDs,sDimsjihm(cid:3)jhƒ(cid:3)ƒe,Me,Minccoohheererenntt(cid:3)(cid:3)ffððDDsismimjjhh,,MMinccoohheererennttÞÞppððhhjjMMcionhceorheenrteÞndtÞhdh: saanmnepalelirngw.aTshseettottaolnthuemsbaemroefvsaamlueplionfgtthriealslawstereesoeftptoop1u0l0a0ti0o0n,
sameasthenumberofparticlesinpopulationannealing.Asshown
Larger BABC indicates that the coherent FFL model is selected
CI inTable 1,themeansofBayesfactorscomputedbyABCrejection
with stronger evidence against the incoherent FFL model. In the
samplerweresimilartothosecomputedbypopulationannealing.
current test study, model selection was conducted with the
Inthecaseoflargeormiddlevaluesofe,thestandarddeviationsof
observed data generated from the coherent FFL model (Fig-
ure 1.C).Thus,BABC shouldbealargevalue,atleastlargerthan Bayes factors computed by ABC rejection sampler were also
CI similar to those of computed by population annealing. However,
1toproposethevalidityofBayesianmodelselectionbypopulation
inthecaseofsmallvalueofe,thestandarddeviationoftheBayes
annealing.
Computation of BABC was done under different annealing factorcomputedbyABCrejectionsamplerwaslargerthanthatof
CI computed by population annealing. Small standard deviation
schedulesinpopulationannealing.Annealingschedulesweresetto
indicates a stable computational result. Thus, these results
same as those in parameter inference. For each annealing
demonstrate the efficiency of population annealing for computa-
schedule,themeanandthestandarddeviationof10independent
computationsofBABC wereshowninTable1.InTable1,asthe tion oftheBayes factors.
CI
PLOSONE | www.plosone.org 10 August2014 | Volume 9 | Issue 8 | e104057