We extend a recently proposed gradient-matching method for inferring interactions in complex systems described by differential equations in various respects: improved gradient inference, evaluation of the influence of the prior on kinetic parameters, comparative evaluation of two model selection paradigms: marginal likelihood versus DIC (divergence information criterion), comparative evaluation of different numerical procedures for computing the marginal likelihood, extension of the methodology from protein phosphorylation to transcriptional regulation, based on a realistic simulation of the underlying molecular processes with Markov jump processes

Aderhold, Andrej

Grzegorczyk, Marco

Husmeier, Dirk

English

Enlighten: Publications

Network reconstruction with realistic models
Marco Grzegorczyk1, Andrej Aderhold2, Dirk Husmeier2
1 Johann Bernoulli Institute, Groningen University, The Netherlands
2 School of Mathematics and Statistics, Glasgow University, UK
E-mail for correspondence: m.a.grzegorczyk@rug.nl
Abstract: We extend a recently proposed gradient-matching method for infer-
ring interactions in complex systems described by differential equations in various
respects: improved gradient inference, evaluation of the influence of the prior on
kinetic parameters, comparative evaluation of two model selection paradigms:
marginal likelihood versus DIC (divergence information criterion), comparative
evaluation of different numerical procedures for computing the marginal likeli-
hood, extension of the methodology from protein phosphorylation to transcrip-
tional regulation, based on a realistic simulation of the underlying molecular
processes with Markov jump processes.
Keywords: Gene Networks; Semi-Mechanistic Models; Bayesian Model Selection
1 INTRODUCTION
A challenging problem for computational statistics is to infer the structure
of regulatory networks from postgenomic data. Two approaches can be dis-
tinguished. The first paradigm aims to apply generic models like Bayesian
networks. The second paradigm is based on mechanistic models and the
detailed mathematical description of the underlying interactions with dif-
ferential equations (DEs). The advantage of this paradigm is a more faith-
ful representation of the interactions. The disadvantage are the substan-
tially higher computational costs of inference. A novel approach, presented
by Oates et al. (2014) and termed “chemical model averaging” (CheMa),
aims for a compromise between the two paradigms by gradient matching.
Given the concentration time series of some quantities (“species”) whose
interactions are to be inferred, the temporal derivatives of the concentra-
tions are estimated from the data. These derivatives are then matched
against those predicted from the DEs by standard statistical techniques.
This paper was published as a part of the proceedings of the 30th Interna-
tional Workshop on Statistical Modelling, Johannes Kepler Universita¨t Linz, 6–10
July 2015. The copyright remains with the author(s). Permission to reproduce or
extract any parts of this abstract should be requested from the author(s).
2 Network reconstruction based on semi-mechanistic models
The resulting “semi-mechanistic” model is effectively a non-linear regres-
sion model, whose computational complexity of inference sits between the
two paradigms. We take the work of Oates et al. (2014) further in five
respects. We improve the gradient matching by adopting techniques from
nonparametric Bayesian statistics with Gaussian processes. We assess the
influence of the parameter prior. We carry out a comparative evaluation
study to assess two model selection paradigms – the marginal likelihood
versus the divergence information criterion. We adapt the method to a new
type of application, namely to transcriptional regulation from gene expres-
sion data. We evaluate the method in a realistic simulation study based on
a stochastic process description of the underlying molecular processes.
2 METHODS
A biopathway can be modelled as a system of ordinary DEs:
dxi
dt
|t=tj = ci − v0,ixi(tj) + fi(xi(tj),θ) (1)
where i is one of N species, xi(tj) is the concentration of i at time point
tj (j = 1, . . . , T ), ci is a baseline production rate, v0,i is a decay rate,
fi(.) is a regulation function, θ are parameters, and xi(tj) is a vector of
concentrations of species that regulate species i. We follow Oates et al.
(2014) and estimate the time derivatives dxidt |t=tj from the observed data
D, and treat the problem as non-linear regression with the likelihood:
p(D|θ) =
N∏
i=1
T∏
j=1
N (ξi(tj)|fi(xi(tj),θ)− v0,ixi(tj), σ2i ) (2)
where ξi(tj) =
dxi
dt |t=tj , and N (.|µ, σ2) is the PDF of a normal distribution.
Oates et al. (2014) obtained the temporal derivatives ξ(tj) with a finite
difference quotient (“numerical gradient”). We here propose to apply a
Gaussian process to smooth interpolation (“analytical gradient”). The
superiority of this approach was recently established by Aderhold et al.
(2014). In the CheMA approach Eq. (1) is implemented with ci = 0, and
fi(.) describes Michaelis-Menten kinetics:
ξi(tj) =
dxi
dt
|t=tj = −v0,ixi(tj) +
∑
j∈pii
vj,i
Ij,i · xj(tj) + (1− Ij,i) · kj,i
xj(tj) + kj,i
(3)
where the sum is over all species j ∈ pii that are regulators of i. The indi-
cator function Ij,i indicates whether species j is an activator or inhibitor.
The term −v0,ixi(t) takes the degradation of xi(t) into account, while the
parameters vj,i and kj,i are the “maximum reaction rate” and “Michaelis-
Menten” parameters. Oates et al. (2004) impose truncated Normal dis-
tributions on kj,i, kj,i ∼ N{kj,i≥0}(1, ν), where ν > 0, use Jeffrey’s prior
Grzegorczyk et al. 3
for σ2i , and a truncated g-prior on the vector Vi of maximum reaction
rate parameters: Vi ∼ N{Vi≥0}(1, Tσ2i (DTiDi)−1), where Di is the design-
matrix for species i. We show that a truncated ridge regression prior:
Vi ∼ N{Vi≥0}(1, δ2i σ2i I), where I is the identity matrix, and δ2i has an
inverse Gamma prior, δ2i ∼ IG(aδ, bδ), yields a better network reconstruc-
tion accuracy. We refer to these methods as CheMa (with the g-prior on
Vi) and iCheMa (improved CheMa with the ridge regression prior on Vi ).
To infer the regulator sets pii of the interaction processes described by
Eq. (1), we compare the divergence information criterion (DIC) and the
marginal likelihood (MLL). DIC is defined as
DIC(pii) = 2 log p(D|θ, pii)− 4
∫
log p(D|θ, pii)p(θ|pii, D)dθ
where θ =
∫
θp(θ|pii, D)dθ is the posterior mean of the parameters. The
integrals are approximated by sums over parameters sampled from the
posterior distribution p(θ|G,D) with MCMC. The marginal likelihood is
p(D|pii) =
∫
p(D|θ, pii)p(θ|pii)dθ (4)
We compare two methods for approximating Eq. (4): Chib’s method and
thermodynamic integration (TI). Chib’s method is based on
p(D|pii) = p(D|θ
?, G)p(θ?|pii)
p(θ?|pii, D) (5)
where the posterior near a selected parameters θ?, p(θ?|G,D), is approxi-
mated with MCMC. TI is based on the power posteriors:
p(θ|pii, D, τ) = p(D|θ, pii)
τp(θ|pii)∫
p(D|θ′, pii)τp(θ′|pii)dθ′
(6)
from which the marginal likelihood is computed via
p(D|pii) =
∫ 1
0
Eθ,τ
[
log p(D|θ, pii)
]
dτ (7)
Here Eθ,τ [.] is an expectation w.r.t. the power posterior in Eq. (6). These
expectations are computed for various temperatures 0 ≤ τ ≤ 1 with popu-
lation MCMC, and the integral in Eq. (7) is then approximated with the
trapezium sum. We choose 10 temperatures τi = (
i
9 )
m (0 ≤ i ≤ 9), and we
vary the exponent m ∈ {4, 8} to obtain TI-4 and TI-8.
3 DATA
We generate T = 240 data points for four species x1, . . . , x4 from iidN (0, 1)
distributions. Subsequently, to obtain non-negative concentrations, the ob-
servations of each individual species are shifted such that the lowest value
4 Network reconstruction based on semi-mechanistic models
0.01 0.1 1 10 100 10000 1e+08 1e+16 1e+20
−30
−20
−10
0
10
20
30
40
50
60
spread factor
D
iff
er
en
ce
 to
 tr
ue
 p
ar
en
t s
et
up
 
 
DIC
Chib
Chib−naiv
TI−4
TI−8
0.01 0.1 1 10 100 10000 1e+08 1e+16 1e+20
−250
−200
−150
−100
−50
0
50
100
150
spread factor
D
iff
er
en
ce
 to
 tr
ue
 p
ar
en
t s
et
up
 
 
DIC
Chib
Chib−naiv
TI−4
TI−8
FIGURE 1. Performance of iCheMa on synthetic data The box plots show
the log marginal likelihood differences between the true regulator set piy = {x2}
and an alternative over-complex set that includes one irrelevant regulator. The
panels correspond to different prior distributions of the parameters. Left panel:
the prior variance of kj,i was kept fixed at ν = 0.5, and the hyperparameter δ
2
i
of the truncated ridge regression prior on Vi was set to spread factor. Right
panel: both hyperparameters δ2i and ν were set to spread factor. The box plots
show the distributions of the average MLL differences across 9 parameter pairs
(v0,y, v2,y), computed for 10 independent data instantiations. Positive differences
indicate that the true model is favoured. Each panel shows: (i) the DIC difference
(DIC), the MLL differences, approximated with (ii) a naive implementation of
Chib’s method (Chib-naiv), (iii) an improved implementation of Chib’s method
(Chib), (iv) TI with m = 4 (TI-4), and (v) TI with m = 8 (TI-8).
is equal to 0, before we re-scale the observations of each species to mean
1. With x1 taking the role of the degradation process and x2 being an
activating regulator (I2,y = 0) of a gradient ξy, which we here assume
to be directly observable, we generate target observations with Eq. (3):
ξy(tj) = −v0,yx1(tj)+v2,y x2(tj)x2(tj)+k2,y + tj , where tj ∼ N (0, σ2) is additive
iid Gaussian noise. We keep k2,y = 1 fixed, and vary the rates (v0,y, v2,y) =
{(1, 1), (0.5, 1), (1.5, 1), (2, 1), (0.2, 1), (2, 0.2), (3, 0.1), (0.2, 2), (0.1, 2)}. Our
goal is to infer the true regulator set piy = {x2} out of all subsets of
{x2, x3, x4}. The degradation, modelled via x1, is included in all models.
We also use the benchmark data from Aderhold et al. (2014), which contain
realistically simulated gene expression time series for genes in the circadian
clock of Arabidopsis thaliana. We focus on the mRNa data points and those
time series generated for the wildtype circadian gene network, shown in the
left panel of Figure 2. The molecular interactions were modelled as individ-
ual discrete events with a Markov jump process and practically simulated
with the Biopepa software, see Aderhold et al. (2014) for details.
Finally, we apply the improved CheMa model to Arabidopsis gene expres-
sion data which were recently measured under the EU-FP7-funded Timing
Metabolism (TiMet) research project. For space restrictions the latter re-
sults are not shown in this paper.
Grzegorczyk et al. 5
0.70 0.72 0.74 0.76 0.78 0.80 0.82
0.
68
0.
70
0.
72
0.
74
AUROC
AU
PR
EC
iCheMA
analytic gradient
iCheMA
numeric gradient
CheMA
analytic gradient
CheMA
numeric gradient
FIGURE 2. Network reconstruction accuracy for the realistic wildtype
Arabdidopsis gene network, based on marginal interaction probabili-
ties (“model averaging”). Left panel: Hypothetical circadian clock network in
A.thaliana from Pokhilko et al. (2010). Right panel: Mean AUROC and AUPREC
scores to quantify the effects of the prior (’CheMa vs. iCheMa’) and the gradi-
ent type (’numerical vs. analytical’) on the network reconstruction accuracy.
4 RESULTS
SYNTHETIC DATA: Figure 1 shows the DIC and log marginal likeli-
hood (MLL) differences between the true and an over-complex regulator
set for different parameter priors. Each parameter prior was a Gaussian
centred on µ = 1, with different variances. For low variances, both DIC
and MLL clearly favour the true network, because the prior ‘pulls’ the spu-
rious interaction parameter from its true value of zero towards a wrong
value of µ = 1. As the prior becomes more diffuse, both the DIC and MLL
differences become less pronounced, but still select the true model up to
spread factors of about 100. As the prior becomes more diffuse, with the
spread factor exceeding 100, DIC fails to select the correct model. MLL, on
the other hand, starts to increasingly favour the true model as the spread
factor further increases beyond 1000. This is a consequence of Lindley’s
paradox, whereby MLL increasingly penalizes the over-complex model for
increasingly vague priors.
The left panel of Figure 1 shows that the different ways of computing
the MLL give very similar results up to a prior spread factor of about
1e+08. For spread factors exceeding this value, the results differ. The MLL
computed with Chib’s method monotonically increases, as expected from
Lindley’s paradox. MLL computed with TI reaches a plateau, with differ-
ent values obtained for different temperature schemes (m = 4 and m = 8).
This is a numerical discretization error that results from the form of the
integrand in Eq. (7), which has most of its area concentrated on values near
τ = 0. The right panel shows that a naive implementation of Chib’s method
can run into numerical instabilities, as diffuse prior distributions can yield
6 Network reconstruction based on semi-mechanistic models
suboptimal attractor states in parameter configuration space. We fixed this
instability by selecting exclusively pivot parameter vectors θ? that are rep-
resentative for the sampled parameter values. (We never observed unstable
results when selecting the parameter vector θ? with the highest posterior
probability from only the sample phase of the MCMC simulation, rather
than the total MCMC trajectory with the burn-in phase included.) With
this numerical stabilization, Chib’s method, whose numerical complexity
compared to TI is lower by a factor of about 10 (as we use K = 10 tem-
peratures τi for TI), is the favourite method for the realistic network data.
REALISTIC NETWORK DATA: The right panel of Figure 2 shows
average AUROC (area under the ROC curve) and AUPREC (area under
the precision recall curve) scores, obtained for five realistically simulated
time series of the wildtype circadian clock network (shown in the left panel).
The proposed analytical gradient (based on a Gaussian process) yields a
significantly improved network reconstruction accuracy for both models:
(i) CheMa with the truncated g-prior on Vi, as proposed by Oates et al.
(2014) and (ii) iCheMa with the truncated ridge regression prior, intro-
duced here. It can also be seen that iCheMa has a slightly better network
prediction accuracy than the original CheMa method given the numerical
gradient, and a substantially better accuracy for the analytic gradient.
5 FURTHER RESULTS
At IWSM 2015 more results will be presented. For space restrictions some
results could not be included in this paper. Most importantly, we will (i)
show that the iCheMa model is superior to established network reconstruc-
tion methods, such as Hierarchical Bayesian regression, Sparse regression
with L1 penalty (Lasso), Sparse regression with L1 and L2 penalty (Elas-
ticNet), Sparse regression with change-points (Tesla), Sparse Bayesian re-
gression, Graphical Gaussian models, Bayesian spline autoregression, State-
space models, Gaussian processes, and Bayesian networks, and we will (ii)
apply the novel iCheMa method to reverse-engineer the circadian clock
network in A. thaliana from TiMet gene expression time series.
References
Aderhold, A., Husmeier, D., and Grzegorczyk, M. (2014). Statistical infer-
ence of regulatory networks for circadian regulation. Statistical appli-
cations in genetics and molecular biology, 13(3), 227 – 273.
Oates et al. (2014) Causal network inference using biochemical kinetics.
Bioinformatics, 30(17), 468 – 474.
Pokhilko et al. (2010). Data assimilation constrains new connections and
components in a complex, eukaryotic circadian clock model. Molecu-
lar Systems Biology, 6(1), online article.


Network Reconstruction with Realistic Models

Enlighten

Network reconstruction with realistic modelsMarco Grzegorczyk1, Andrej Aderhold2, Dirk Husmeier21 Johann Bernoulli Institute, Groningen University, The Netherlands2 School of Mathematics and Statistics, Glasgow University, UKE-mail for correspondence: m.a.grzegorczyk@rug.nlAbstract: We extend a recently proposed gradient-matching method for infer-ring interactions in complex systems described by differential equations in variousrespects: improved gradient inference, evaluation of the influence of the prior onkinetic parameters, comparative evaluation of two model selection paradigms:marginal likelihood versus DIC (divergence information criterion), comparativeevaluation of different numerical procedures for computing the marginal likeli-hood, extension of the methodology from protein phosphorylation to transcrip-tional regulation, based on a realistic simulation of the underlying molecularprocesses with Markov jump processes.Keywords: Gene Networks; Semi-Mechanistic Models; Bayesian Model Selection1 INTRODUCTIONA challenging problem for computational statistics is to infer the structureof regulatory networks from postgenomic data. Two approaches can be dis-tinguished. The first paradigm aims to apply generic models like Bayesiannetworks. The second paradigm is based on mechanistic models and thedetailed mathematical description of the underlying interactions with dif-ferential equations (DEs). The advantage of this paradigm is a more faith-ful representation of the interactions. The disadvantage are the substan-tially higher computational costs of inference. A novel approach, presentedby Oates et al. (2014) and termed “chemical model averaging” (CheMa),aims for a compromise between the two paradigms by gradient matching.Given the concentration time series of some quantities (“species”) whoseinteractions are to be inferred, the temporal derivatives of the concentra-tions are estimated from the data. These derivatives are then matchedagainst those predicted from the DEs by standard statistical techniques.This paper was published as a part of the proceedings of the 30th Interna-tional Workshop on Statistical Modelling, Johannes Kepler Universita¨t Linz, 6–10July 2015. The copyright remains with the author(s). Permission to reproduce orextract any parts of this abstract should be requested from the author(s).2 Network reconstruction based on semi-mechanistic modelsThe resulting “semi-mechanistic” model is effectively a non-linear regres-sion model, whose computational complexity of inference sits between thetwo paradigms. We take the work of Oates et al. (2014) further in fiverespects. We improve the gradient matching by adopting techniques fromnonparametric Bayesian statistics with Gaussian processes. We assess theinfluence of the parameter prior. We carry out a comparative evaluationstudy to assess two model selection paradigms – the marginal likelihoodversus the divergence information criterion. We adapt the method to a newtype of application, namely to transcriptional regulation from gene expres-sion data. We evaluate the method in a realistic simulation study based ona stochastic process description of the underlying molecular processes.2 METHODSA biopathway can be modelled as a system of ordinary DEs:dxidt|t=tj = ci − v0,ixi(tj) + fi(xi(tj),θ) (1)where i is one of N species, xi(tj) is the concentration of i at time pointtj (j = 1, . . . , T ), ci is a baseline production rate, v0,i is a decay rate,fi(.) is a regulation function, θ are parameters, and xi(tj) is a vector ofconcentrations of species that regulate species i. We follow Oates et al.(2014) and estimate the time derivatives dxidt |t=tj from the observed dataD, and treat the problem as non-linear regression with the likelihood:p(D|θ) =N∏i=1T∏j=1N (ξi(tj)|fi(xi(tj),θ)− v0,ixi(tj), σ2i ) (2)where ξi(tj) =dxidt |t=tj , and N (.|µ, σ2) is the PDF of a normal distribution.Oates et al. (2014) obtained the temporal derivatives ξ(tj) with a finitedifference quotient (“numerical gradient”). We here propose to apply aGaussian process to smooth interpolation (“analytical gradient”). Thesuperiority of this approach was recently established by Aderhold et al.(2014). In the CheMA approach Eq. (1) is implemented with ci = 0, andfi(.) describes Michaelis-Menten kinetics:ξi(tj) =dxidt|t=tj = −v0,ixi(tj) +∑j∈piivj,iIj,i · xj(tj) + (1− Ij,i) · kj,ixj(tj) + kj,i(3)where the sum is over all species j ∈ pii that are regulators of i. The indi-cator function Ij,i indicates whether species j is an activator or inhibitor.The term −v0,ixi(t) takes the degradation of xi(t) into account, while theparameters vj,i and kj,i are the “maximum reaction rate” and “Michaelis-Menten” parameters. Oates et al. (2004) impose truncated Normal dis-tributions on kj,i, kj,i ∼ N{kj,i≥0}(1, ν), where ν > 0, use Jeffrey’s priorGrzegorczyk et al. 3for σ2i , and a truncated g-prior on the vector Vi of maximum reactionrate parameters: Vi ∼ N{Vi≥0}(1, Tσ2i (DTiDi)−1), where Di is the design-matrix for species i. We show that a truncated ridge regression prior:Vi ∼ N{Vi≥0}(1, δ2i σ2i I), where I is the identity matrix, and δ2i has aninverse Gamma prior, δ2i ∼ IG(aδ, bδ), yields a better network reconstruc-tion accuracy. We refer to these methods as CheMa (with the g-prior onVi) and iCheMa (improved CheMa with the ridge regression prior on Vi ).To infer the regulator sets pii of the interaction processes described byEq. (1), we compare the divergence information criterion (DIC) and themarginal likelihood (MLL). DIC is defined asDIC(pii) = 2 log p(D|θ, pii)− 4∫log p(D|θ, pii)p(θ|pii, D)dθwhere θ =∫θp(θ|pii, D)dθ is the posterior mean of the parameters. Theintegrals are approximated by sums over parameters sampled from theposterior distribution p(θ|G,D) with MCMC. The marginal likelihood isp(D|pii) =∫p(D|θ, pii)p(θ|pii)dθ (4)We compare two methods for approximating Eq. (4): Chib’s method andthermodynamic integration (TI). Chib’s method is based onp(D|pii) = p(D|θ?, G)p(θ?|pii)p(θ?|pii, D) (5)where the posterior near a selected parameters θ?, p(θ?|G,D), is approxi-mated with MCMC. TI is based on the power posteriors:p(θ|pii, D, τ) = p(D|θ, pii)τp(θ|pii)∫p(D|θ′, pii)τp(θ′|pii)dθ′(6)from which the marginal likelihood is computed viap(D|pii) =∫ 10Eθ,τ[log p(D|θ, pii)]dτ (7)Here Eθ,τ [.] is an expectation w.r.t. the power posterior in Eq. (6). Theseexpectations are computed for various temperatures 0 ≤ τ ≤ 1 with popu-lation MCMC, and the integral in Eq. (7) is then approximated with thetrapezium sum. We choose 10 temperatures τi = (i9 )m (0 ≤ i ≤ 9), and wevary the exponent m ∈ {4, 8} to obtain TI-4 and TI-8.3 DATAWe generate T = 240 data points for four species x1, . . . , x4 from iidN (0, 1)distributions. Subsequently, to obtain non-negative concentrations, the ob-servations of each individual species are shifted such that the lowest value4 Network reconstruction based on semi-mechanistic models0.01 0.1 1 10 100 10000 1e+08 1e+16 1e+20−30−20−100102030405060spread factorDifference to true parent setup  DICChibChib−naivTI−4TI−80.01 0.1 1 10 100 10000 1e+08 1e+16 1e+20−250−200−150−100−50050100150spread factorDifference to true parent setup  DICChibChib−naivTI−4TI−8FIGURE 1. Performance of iCheMa on synthetic data The box plots showthe log marginal likelihood differences between the true regulator set piy = {x2}and an alternative over-complex set that includes one irrelevant regulator. Thepanels correspond to different prior distributions of the parameters. Left panel:the prior variance of kj,i was kept fixed at ν = 0.5, and the hyperparameter δ2iof the truncated ridge regression prior on Vi was set to spread factor. Rightpanel: both hyperparameters δ2i and ν were set to spread factor. The box plotsshow the distributions of the average MLL differences across 9 parameter pairs(v0,y, v2,y), computed for 10 independent data instantiations. Positive differencesindicate that the true model is favoured. Each panel shows: (i) the DIC difference(DIC), the MLL differences, approximated with (ii) a naive implementation ofChib’s method (Chib-naiv), (iii) an improved implementation of Chib’s method(Chib), (iv) TI with m = 4 (TI-4), and (v) TI with m = 8 (TI-8).is equal to 0, before we re-scale the observations of each species to mean1. With x1 taking the role of the degradation process and x2 being anactivating regulator (I2,y = 0) of a gradient ξy, which we here assumeto be directly observable, we generate target observations with Eq. (3):ξy(tj) = −v0,yx1(tj)+v2,y x2(tj)x2(tj)+k2,y + tj , where tj ∼ N (0, σ2) is additiveiid Gaussian noise. We keep k2,y = 1 fixed, and vary the rates (v0,y, v2,y) ={(1, 1), (0.5, 1), (1.5, 1), (2, 1), (0.2, 1), (2, 0.2), (3, 0.1), (0.2, 2), (0.1, 2)}. Ourgoal is to infer the true regulator set piy = {x2} out of all subsets of{x2, x3, x4}. The degradation, modelled via x1, is included in all models.We also use the benchmark data from Aderhold et al. (2014), which containrealistically simulated gene expression time series for genes in the circadianclock of Arabidopsis thaliana. We focus on the mRNa data points and thosetime series generated for the wildtype circadian gene network, shown in theleft panel of Figure 2. The molecular interactions were modelled as individ-ual discrete events with a Markov jump process and practically simulatedwith the Biopepa software, see Aderhold et al. (2014) for details.Finally, we apply the improved CheMa model to Arabidopsis gene expres-sion data which were recently measured under the EU-FP7-funded TimingMetabolism (TiMet) research project. For space restrictions the latter re-sults are not shown in this paper.Grzegorczyk et al. 50.70 0.72 0.74 0.76 0.78 0.80 0.820.680.700.720.74AUROCAUPRECiCheMAanalytic gradientiCheMAnumeric gradientCheMAanalytic gradientCheMAnumeric gradientFIGURE 2. Network reconstruction accuracy for the realistic wildtypeArabdidopsis gene network, based on marginal interaction probabili-ties (“model averaging”). Left panel: Hypothetical circadian clock network inA.thaliana from Pokhilko et al. (2010). Right panel: Mean AUROC and AUPRECscores to quantify the effects of the prior (’CheMa vs. iCheMa’) and the gradi-ent type (’numerical vs. analytical’) on the network reconstruction accuracy.4 RESULTSSYNTHETIC DATA: Figure 1 shows the DIC and log marginal likeli-hood (MLL) differences between the true and an over-complex regulatorset for different parameter priors. Each parameter prior was a Gaussiancentred on µ = 1, with different variances. For low variances, both DICand MLL clearly favour the true network, because the prior ‘pulls’ the spu-rious interaction parameter from its true value of zero towards a wrongvalue of µ = 1. As the prior becomes more diffuse, both the DIC and MLLdifferences become less pronounced, but still select the true model up tospread factors of about 100. As the prior becomes more diffuse, with thespread factor exceeding 100, DIC fails to select the correct model. MLL, onthe other hand, starts to increasingly favour the true model as the spreadfactor further increases beyond 1000. This is a consequence of Lindley’sparadox, whereby MLL increasingly penalizes the over-complex model forincreasingly vague priors.The left panel of Figure 1 shows that the different ways of computingthe MLL give very similar results up to a prior spread factor of about1e+08. For spread factors exceeding this value, the results differ. The MLLcomputed with Chib’s method monotonically increases, as expected fromLindley’s paradox. MLL computed with TI reaches a plateau, with differ-ent values obtained for different temperature schemes (m = 4 and m = 8).This is a numerical discretization error that results from the form of theintegrand in Eq. (7), which has most of its area concentrated on values nearτ = 0. The right panel shows that a naive implementation of Chib’s methodcan run into numerical instabilities, as diffuse prior distributions can yield6 Network reconstruction based on semi-mechanistic modelssuboptimal attractor states in parameter configuration space. We fixed thisinstability by selecting exclusively pivot parameter vectors θ? that are rep-resentative for the sampled parameter values. (We never observed unstableresults when selecting the parameter vector θ? with the highest posteriorprobability from only the sample phase of the MCMC simulation, ratherthan the total MCMC trajectory with the burn-in phase included.) Withthis numerical stabilization, Chib’s method, whose numerical complexitycompared to TI is lower by a factor of about 10 (as we use K = 10 tem-peratures τi for TI), is the favourite method for the realistic network data.REALISTIC NETWORK DATA: The right panel of Figure 2 showsaverage AUROC (area under the ROC curve) and AUPREC (area underthe precision recall curve) scores, obtained for five realistically simulatedtime series of the wildtype circadian clock network (shown in the left panel).The proposed analytical gradient (based on a Gaussian process) yields asignificantly improved network reconstruction accuracy for both models:(i) CheMa with the truncated g-prior on Vi, as proposed by Oates et al.(2014) and (ii) iCheMa with the truncated ridge regression prior, intro-duced here. It can also be seen that iCheMa has a slightly better networkprediction accuracy than the original CheMa method given the numericalgradient, and a substantially better accuracy for the analytic gradient.5 FURTHER RESULTSAt IWSM 2015 more results will be presented. For space restrictions someresults could not be included in this paper. Most importantly, we will (i)show that the iCheMa model is superior to established network reconstruc-tion methods, such as Hierarchical Bayesian regression, Sparse regressionwith L1 penalty (Lasso), Sparse regression with L1 and L2 penalty (Elas-ticNet), Sparse regression with change-points (Tesla), Sparse Bayesian re-gression, Graphical Gaussian models, Bayesian spline autoregression, State-space models, Gaussian processes, and Bayesian networks, and we will (ii)apply the novel iCheMa method to reverse-engineer the circadian clocknetwork in A. thaliana from TiMet gene expression time series.ReferencesAderhold, A., Husmeier, D., and Grzegorczyk, M. (2014). Statistical infer-ence of regulatory networks for circadian regulation. Statistical appli-cations in genetics and molecular biology, 13(3), 227 – 273.Oates et al. (2014) Causal network inference using biochemical kinetics.Bioinformatics, 30(17), 468 – 474.Pokhilko et al. (2010). Data assimilation constrains new connections andcomponents in a complex, eukaryotic circadian clock model. Molecu-lar Systems Biology, 6(1), online article.

http://eprints.gla.ac.uk/108338/1/108338.pdf

Network Reconstruction with Realistic Models

Abstract

Similar works

Full text

Available Versions

Enlighten: Publications

Enlighten