BACKGROUND
the availability of high throughput experimental platforms has transformed molecular biology into a data-rich science.
however, the development of computational and mathematical tools to extract and interpret the wealth of information hidden in these data is lagging behind.
for example, genome wide rna interference screens have enabled the phenotypic characterization of genes at an unprecedented scale in living cells
<cit> .
however, the interpretation of such data and the placement of hits in their functional and temporal context in cellular pathways remains a major challenge
<cit> .
machine learning can address many of the questions arising from the interpretation of large scale molecular biological data, and has become an important tool of bioinformatics and systems biology.
automatic network inference deals with the reconstruction of regulatory or signaling networks directly from experimental data using statistical or machine learning approaches, a field that has received significant attention in the last decade.
there is a wide variety of different approaches available that can be used to infer genetic regulatory or signal transduction networks from experimental data.
approaches employed include bayesian networks
and dynamic bayesian networks
, boolean models
, auto-regressive models, correlation-based and mutual-information based models, clustering techniques, differential equation models, and others
.
these methods differ in the level of detail at which they reconstruct networks, and in their underlying assumptions and data requirements.
some produce an undirected graph, where edges do not indicate which gene or protein in a connected pair is the activator, and which is the activated gene or protein
<cit> .
others specify the regulator with directed edges
<cit> , and a few label the edges with kinetic parameters
<cit> .
there are a number of successful applications of network inference approaches to elucidate cellular signaling pathways, including meta-approaches integrating different methods
<cit> .
among the successful applications, for example, sachs et al.
used bayesian networks to reconstruct cellular protein signaling networks from protein phosphorylation measurements
<cit> .
nelander et al.
used a model based on nonlinear differential equations to infer signaling networks in cancer from combinatorial drug perturbation data
<cit> .
eduati et al.
demonstrate the integration of literature-knowledge into data driven approaches, and show a successful application to a signaling network related to growth-signaling and inflammation
<cit> .
ciaccio et al.
used bayesian networks and two different mutual-information based approaches to infer signaling networks downstream of the egf receptor
<cit> .
hill et al.
use dynamic bayesian networks to study signaling in a cancer cell line
<cit> .
other approaches include nested effects models
<cit> , deterministic effects propagation networks
<cit>  or probabilistic boolean threshold networks
<cit> , and have been applied, for example, to reconstruct signaling networks in the erbb pathway, or in the innate immune response to infection.
network reconstruction can be performed from observational data alone.
however, the quality of the reconstruction increases substantially if experimental perturbations followed by observations of the system’s dynamic response are available
<cit> .
surprisingly, while there are many approaches available for time course data, and several approaches for perturbation data, there are only few methods available that can handle both time course data and perturbations at the same time.
exceptions are, for example, the differential equation approach presented by nelander et al.
<cit> , dynamic nested effects models
<cit>  or dynamic deterministic effects propagation models
<cit> .
among the stochastic approaches, dnems rely on high-dimensional, indirect readouts of rather qualitative knockdown "effects", such as microarrays performed at different time points after every knockdown, or multidimensional features derived from live cell imaging.
such data are often not available, and also provide only indirect information about the underlying signaling pathway.
fröhlich et al.
proposed deterministic effects propagation networks  to reconstruct networks from perturbation data with direct observation of involved proteins, measured e.g.
using reverse phase protein arrays
<cit> .
however, depns treat each time point as an independent measurement and do not model the time dependent behavior of the system explicitly.
bender et al.
proposed dynamic depns , to take the availability of longitudinal data as well as inhibitory interactions in a network into account
<cit> .
d-depns use the viterbi algorithm to identify the state transitions of a hidden markov model from the time course data, with states corresponding to combinations of activities of nodes in the network.
a likelihood function is defined to score network models given the estimated state transitions, and the network space is then searched using a genetic algorithm.
due to this underlying procedure, d-depns require relatively long time series , and have substantial running time requirements if large numbers of different perturbations involving most or all of the proteins in a network are performed.
we here focus on the problem of reconstructing signaling networks from short time series with only two or three time points, after a large number of different, possibly combinatorial perturbations, targeting most or all of the genes in the network.
such data arise, for example, if rnai experiments are combined with protein array measurements at few time points after the perturbation, or in time-resolved mass spectrometric assays under different conditions.
we propose a new method for the identification of kinetic models of signaling networks from such data, dynamic probabilistic boolean threshold networks , which can treat multiple, combinatorial interventions as well as incomplete observations.
through the use of a fairly simple, discrete-state model where proteins are either active or inactive, our method is applicable to qualitative data and does not require detailed quantitative measurements.
furthermore, through the bayesian framework employed, we can easily handle noisy or missing data, and using mcmc sampling, can analyze full posterior distributions over model parameters, yielding information about different, alternative network topologies that are consistent with the experimental data.
our work builds upon probabilistic boolean threshold networks  regarding the underlying dynamic network model
<cit> , but extends the method in the following aspects:  we fully exploit the information from time course data, by explicitely integrating time into the model likelihood.
in contrast, pbtns in their original implementation only exploit observational data taken at a single time point, usually at steady state.
pbtns are limited to relatively small networks of around  <dig> nodes.
through a more efficient sampling method, the approach presented here can be used for larger networks, effectively more than doubling the feasible network size.
pbtns were implemented with a focus on data with only single downstream observations at a single time point, and therefore similarly to other  bayesian network based approaches, do not allow feedback loops in the network.
the approach presented here can handle feedback loops in the network.
d-pbtns can handle overexpression and knockdown data, in contrast to pbtns, which are limited to knockdown data.
our approach uses discrete states of the proteins in the network, which can be either "on" or "off".
in terms of the biological interpretation, this could be phosphorylation networks, in which case edges would correspond to phosphorylation or dephosphorylation events.
alternatively, assuming that proteins are "present" or "absent", edges in the model could also be interpreted as transcriptional regulation.
correspondingly, nodes can represent "activated" protein levels or total protein abundance depending on context.
we evaluate the performance of our method on simulated data, and study its behavior on networks of different sizes, with different amounts of available data and different levels of noise, and compare results against pbtn
<cit> , depn
<cit>  and bdagl, a dynamic bayesian network based approach
<cit> .
dynamic depns could not be used for the comparison, as they were unable to handle the short time courses and large number of perturbations.
we show inference results for erbb signaling in breast cancer using depn, pbtn and d-pbtn, demonstrating superior performance of d-pbtn, and showing cross-talk between different branches of the erbb pathway.
methods
mathematical model
we model a signaling network by a weighted, directed graph
.
cycles are permitted.
the node set
corresponds to proteins; an edge
from node vi to vj represents a regulation of vj by vi, such as activation or deactivation via phosphorylation.
the strength and type of an interaction is specified by
, with wi,j >  <dig> for activations, wi,j <  <dig> for inactivations and wi,j :=  <dig> if
.
given n, the number of nodes in the network, we can thus describe the network topology by its n × n weighted adjacency matrix
.
we furthermore consider each node vi a boolean variable.
model time is assumed discrete, and all nodes are updated simultaneously between two time steps.
the state of node vi at model time t +  <dig> is described by a stochastic function of the states of all nodes
at model time t:
<dig>
the relationship between a target protein and its regulators is hence modeled by a sigmoid function, an ansatz for nonlinear systems with saturation phenomena.
equation  activates or deactivates a target protein depending on the weighted sum of incoming regulations.
the level of stochasticity in this process is defined by the positive parameter γ, whereas w <dig> i defines a basal activation probability of vi in the absence of any regulatory effects.
in case of interventions, the state of a node v targeted by perturbation is no longer subject to the stochastic dynamics of the system, instead, its value is fixed through the intervention.
we write
, k ∈ { <dig> …,k}, for the set of nodes targeted simultaneously in perturbation experiment k, and assume that a knockdown of
amounts to fixing all nodes
to the "off" state at all model time points after the knockdown.
similarly, overexpression experiments can be simulated by fixing the affected nodes vi to "on".
let us now assume that we have observed experimentally the states of the nodes vi at t different  time points
,
.
we note in passing that this normally requires a discretization step, as experimental data are typically measured continuously – see the section on the erbb data set below.
let us for now assume that the set of experimental observations
consists of measurements of n discrete protein states at t time points after k different perturbations.
we write
for the observed activity of node vi at the real time point
after knockdown of node set
.
to map real time
to model time t, we introduce a parameter vector
to denote the number of model time steps that correspond to the real time elapsed between the different experimental time points
, that is,
denotes the number of model time steps according to equation  required to transition from real time point
to time point
.
assuming a one-to-one correspondence between model time and real time, that is, assuming that
for all
, we can define a likelihood by calculating the product of
according to equation  over all nodes v, time points t and perturbations k:
<dig>
if τt >  <dig>  one or more intermediate model time steps must be made to transition from experimental time point
to experimental time point
.
as no experimental observations are available for the intermediate time steps, we must marginalize over unobserved time steps in the evaluation of the likelihood .
writing
for the full sequence of model steps, where the sequence of intermediate steps
are unobserved, equation  becomes
<dig>
where
comprises the product over all individual nodes as in equation , and where the sum is over all possible combinations of values for the intermediate model steps.
missing values in the experimental data can be treated similarly – the likelihood can then be computed by integrating equation  over all possible combinations of values for unobserved
.
we note that pbtns are based on a closely related likelihood function, and employ the same underlying model for the relation between nodes, given in equation .
however, pbtns differ in a key aspect from the framework presented here: pbtns do not allow multiple time points in their likelihood, but rather assume that measurements are made at a single time point after knockdown, usually taken at steady state.
inclusion of prior knowledge
in many biological settings, available data are insufficient to unambiguously reconstruct the underlying network.
in such situations, strict regularization of the objective function using for example maximum parsimony
<cit> , or the inclusion of additional prior biological knowledge
<cit>  can help.
both can be done via bayes’ theorem using prior distributions on the model parameters Ω. the model posterior is
<dig>
where
is the likelihood function as defined above, p is an adequately defined prior distribution, and
is a constant normalizing factor that does not depend on Ω and can be neglected in maximization or sampling.
the prior can be written as
, assuming that the different types of model parameters are independent.
then,
describes our belief about the correct topology, prior to seeing any data, and p describes our belief about the speed of signal transduction in the network.
assuming independence of parameters, the full prior can be written as:
<dig>
in the network inference setting, it is unlikely that the true underlying biological network is densely connected.
we rather expect a sparse network, where most of the nodes have only few other nodes acting on them
.
this can be mathematically encoded in a prior on wi,j  <dig>
with positive shape and rate parameters q and s
<cit> .
note that, after taking the negative logarithm and dropping the normalizing factor 1/, this is equivalent to regularization using the lq-norm, and corresponds to a laplacian prior if q =  <dig> or a normal prior if q =  <dig>  and the parameter s determines the width of the distribution.
the parameters μi,j can be used to encode additional knowledge on specific edges, by setting individual μi,j to nonzero values.
if no prior knowledge is available, μi,j defaults to  <dig>  we note that the μi,j can be chosen continuously depending on the expected strength of the effect and certainty in its presence.
we expect signaling molecules to be "off" in the absence of any signal, and therefore use independent negative gamma priors for p,
<dig>
with rate and shape parameters r <dig> and s <dig>  respectively.
we hence assume w <dig> i to be non-positive.
this assumption can easily be replaced if signaling molecules should also be allowed to be "on" in the absence of signal, for example by using a gamma prior on the absolute value |w <dig> i|, or a normal prior on w <dig> i. lastly, as model time is discrete and necessarily positive, we use a binomial prior p on τi,
<dig>
with parameters n and p.
network inference
we could now maximize
to find the most probable network Ω underlying the data
.
this is reasonable if the posterior is unimodal and has a sharp, narrow optimum.
however, since the optimization problem encountered for real data is typically underdetermined, several alternative models Ω often explain the data equally well, and marginal distributions for parameters show wide peaks and large confidence intervals.
we therefore use markov chain monte carlo  sampling to fully characterize the posterior.
the main difficulty in doing this lies in the calculation of the likelihood
, which due to the required marginalization over unobserved nodes and time points quickly becomes very involved and time-consuming, compare equation .
we therefore approximate
by simulating the dynamics of the model with parameters Ω for each knockdown
, and then compare the simulation results mk at each time point t with the experimental data dk.
if this is repeated a large number of times, we can approximate
by the number of times we get
back in the simulation runs.
by combining this with mcmc sampling, a very efficient approach for evaluation of the posterior distribution arises, that does not require an explicit computation of the likelihood
<cit> .
in contrast to pbtn, which employed a metropolis-hastings based sampler, we have now implemented a far more efficient sampler in d-pbtn.
we employ distributed evolutionary monte carlo
<cit> , a sampler that combines features of a distributed genetic algorithm with mcmc sampling.
in brief, demc starts with a population of m markov chains, which are grouped into g subpopulations.
each individual chain describes a specific network Ω. updates within each subpopulation are done using the genetic operators mutation and cross-over, and migration allows individual solutions to move between subpopulations.
thereafter, each chain is scored by initializing each node in the network with the values experimentally observed at time  <dig>  and then simulating using equation .
this is followed by a metropolis-hastings step to accept or reject the new state, to ensure ergodicity of the chains.
a detailed description of the procedure is given in additional file
<dig>
data simulation
since for real experimental data typically no "gold standard" network is available to assess results, we used simulated data to evaluate d-pbtn.
this allows it to systematically alter properties of the data, such as the inherent level of noise or amount of missing values, and to directly assess the performance impact this has.
we simulated data in three different ways:
simulated network  <dig>  is a 7-node feedforward network, shown in figure
1a .
weights wi,j were set to  <dig> for edges shown in the figure, and  <dig> otherwise.
baseline weights were set to w <dig> i = - <dig>  to have a high probability for an unregulated protein to be in the off state.
data simulation for this network was performed assuming that all proteins are in the off state initially, except for the receptor  <dig>  stochastic signaling is then simulated using equation , updating all nodes simultaneously, and using γ =  <dig>  knockdowns were simulated by fixing targeted proteins in the off state.
we simulated knockdowns of all individual proteins, and two combinatorial knockdowns of proteins  and .
two time points were used for network inference, the first one immediately after knockdown and activation of the receptor; the second time point after  <dig> simulation steps, allowing for sufficient time for the signal to propagate through the whole network.figure  <dig>
evaluation on simulated data.
a: the panel shows network topologies used to simulate data.
simulated network  <dig>  is a simple feedforward network, whereas simulated network  <dig>  contains a negative feedback loop.
b: network reconstruction was performed for sn <dig>  changing the number of unobserved proteins.
shown is the distribution of the area under the roc curve  of  <dig> replicates of simulated data sets, over the fraction of unobserved proteins.
from left to right: d-pbtn, depn, bdagl, pbtn.
the dashed line at auc
roc =  <dig>  shows expected results for random guessing.
c: the panel shows the distribution of auc
roc values obtained for different levels of noise on sn <dig>  noise is introduced by switching the state of the indicated fraction of proteins in the "measured" data, thus introducing errors in the data.
d: this panel shows performance of d-pbtn, depn, bdagl and pbtn on networks sampled from kegg, with data simulated as described in methods.
shown are auc
roc values of <dig> simulated data sets, generated from ten different kegg subnetworks.
using the same network topology and parameters, we furthermore simulated overexpression experiments .
overexpressed nodes were fixed to the on state, while we assumed all other nodes to be inactive initially.
we simulated overexpression experiments of all individual nodes, as well as combinatorial overexpression of nodes  and nodes , again by running the simulation over  <dig> time steps as above.
simulated network  <dig>  extends sn <dig> by a feedback loop, see figure
1a .
data for sn <dig> was simulated using the model proposed by fröhlich et al.
<cit> .
in contrast to sn <dig>  signaling is deterministic, and noise arises only at the measurement stage.
we initialized all nodes as above and assumed that a node deterministically becomes active if there are more activating than inhibiting proteins among its parents.
measurements were simulated by sampling from a normal distribution with
-distributed mean for active nodes, and
-distributed mean if the node was inactive.
the variances of the normal distributions were drawn from an inverse scaled chi-squared distribution
, as suggested by fröhlich et al.
<cit>  based on their experimental observations.
measurements were simulated immediately after the knockdown and activation, and after four and six simulation steps.
simulated network  <dig>  comprises a set of ten different networks randomly extracted from the kegg database, each containing seven connected proteins, taking only protein-protein interactions into account.
we simulated time-course data for  <dig> time points as described for sn <dig> for each of the ten subnetworks.
for each subnetwork, in addition to the starting point, two further randomly selected time points were included into the final data set used for network reconstruction.
reported performance results are average performances over all ten subnetworks, avoiding biasing of results towards a single topology or time interval.
implementation and performance evaluation
we implemented d-pbtn in c♯, and evaluated the approach on both simulated and real experimental data.
source code  and additional information is available from our website at
http://www.kaderali.org/software.html.
convergence of the markov chains was assessed using the methods proposed by brooks and gelman
<cit> , results are given in additional file
<dig>  sampled networks were aggregated using halite clustering
<cit> .
halite is a density-based clustering method, that we use to identify clusters of similar networks from the markov chains.
cluster representative networks were computed using the within-cluster median of wi,j.
to quantify network reconstruction performance, we performed receiver operator characteristic  and precision recall  analysis.
in brief, a threshold δ is used for discretization of edge values, where the median
of sampled values for a particular edge wi,j is compared with δ. if
, the edge ei,j is assumed present in the network,
, and absent otherwise,
.
the resulting network with edge set
is then compared against the gold standard network, and sensitivity, specificity and precision are computed for given δ. this is then repeated by varying δ, and sensitivity is plotted over specificity for different δ in a roc plot, and precision over sensitivity in a pr plot.
finally, both roc and pr curves can be summarized by computing the area under the curve  as single numbers summarizing performance over a wide range of threshold values δ. the advantage of this approach is that it is not necessary to pick a specific threshold δ for the discretization, which may be difficult to do as this depends on  user preferences regarding the tradeoff between false positive and false negative edges, but the analysis summarizes performance over the full range of feasible thresholds δ.
we first assessed performance of d-pbtn on simulated data, varying the amount of noise in the data, using different amounts of simulated data and different models for data simulation, and evaluated stability of results for changing model hyperparameters.
we then assessed performance on publicly available data regarding signaling in the erbb network.
inferred networks on real data were assessed using string  <dig>  as reference
<cit> , using only edges with at least 70% confidence in string.
since interactions given in string are undirected and unsigned, we disregarded directional information and edge sign.
we compare performance of our approach with depn
<cit> , the non-time-series version of our approach
<cit> , the bayesian directed acyclic graph learning tool  developed by eaton and murphy
<cit> , and with random guessing.
notably, d-depn could not be used on the simulated data due to the nature of the short time courses used.
results for pbtn were obtained using a modified version of the original pbtn implementation
<cit>  that optimizes the posterior distribution instead of sampling, thus allowing for a more efficient computation on networks with a unimodal posterior.
this compromise is necessary, as the original sampling-version of pbtn is too slow to allow a rigorous evaluation on networks of the size used here.
time points were used as independent measurements in pbtn.
to obtain results for the depns, we used greedy hill-climbing and bootstrapping with  <dig> bootstrap samples, as proposed by the authors of the method.
the bootstrap samples were used to obtain weights on the edges, which were subjected to roc and pr analysis.
we note that depns operate on equivalence classes of graphs, and can not distinguish between graphs within an equivalence class.
we therefore used the transitively closed network as a representative of the full equivalence class to evaluate the depn results.
this is in contrast to the other methods, which in principle are able to determine a unique network, provided sufficient data is given, and which we compared directly to the gold standard network.
lastly, results for bdagl were obtained using the modified marginal likelihood method for perfect interventional data, with uniform prior.
RESULTS
as a first test of our method, we assessed the performance of d-pbtn on simulated data, without noise and with full observations of all nodes, using the network topology sn <dig>  to compare performance of d-pbtn with other approaches, we used the non-time series version pbtn, depn and bdagl.
for d-pbtn, shape and rate parameters of the prior distribution  on w were set to q,s =  <dig>  respectively, corresponding to a standard l <dig> regularization of the edge weights as the simplest conceivable setting.
rate and shape parameters of the negative gamma prior on w <dig> i were set to r0 =  <dig>  s0 =  <dig>  resulting in a moderate deactivation of unregulated proteins, and a binomial prior with n =  <dig> p =  <dig>  was used for τ, giving an expected value of τ of  <dig> time steps.
the stochasticity parameter γ was set to  <dig> , corresponding to an average level of noise in the data.
three sub-populations with four markov chains were used with  <dig> steps each and a burn-in of  <dig>  steps.
parameters for pbtn were set accordingly, using the two time points as independent replicate measurements.
results for all four approaches are shown in table
<dig>  values shown are median values of the aucroc and aucpr measures over  <dig> replicate simulations.
with no noise and no missing values, d-pbtn shows superior performance over the other three approaches on the simulated data.
notably, the comparison between d-pbtn and pbtn shows the added value of the temporal information, indicating that while a significant part of the information in the data comes from the knockdowns, even short time courses with just  <dig> time points are of added value.
bdagl shows weaker results than the other three approaches, both in terms of the aucroc and aucpr, possibly due to a significantly different underlying mathematical model.table  <dig>
performance comparison of d-pbtn, pbtn, depn and bdagl on simulated data
auc
roc
auc
pr
network
sn1
sn2
sn3
shown are achieved values of the area under the roc curve  and the area under the pr curve .
values shown are median values over  <dig> iterations.
inference was performed on the sn <dig> data, with data with different levels of noise and missing values, on the sn <dig> data including a negative feedback loop, and on the sn <dig>  data set.
note that bdagl cannot handle missing values.
the upper part of the table shows the area under the roc curve, while the lower part of the table shows the area under the pr curve.
effect of missing data
we next assessed the effect of missing observations on network inference, using the sn <dig> data.
results are shown in figure
1b for aucroc and additional file
2: figure s <dig> for aucpr, and are summarized in table
<dig>  note that bdagl could not be used in this comparison, as the method cannot handle missing values.
we randomly selected 16%, 33% and 50% of the nodes, and completely removed the corresponding simulated observations for these nodes before carrying out the network inference.
model parameters were set as described above.
computation took approximately  <dig>  hours for d-pbtn on a  <dig> ghz intel 64-bit processor .
for comparison, runtime for depn was only around  <dig> minutes.
regarding quality of the inferred networks, d-pbtn shows superior performance over depn, both with respect to aucroc and aucpr.
pbtn did not perform better than random guessing on all runs with missing values .
both d-pbtn and depn show a decrease in the quality of the inferred network with increasing amounts of missing observations.
in spite of this, even with 50% of missing data, both approaches still perform significantly better than guessing .
median aucroc decreased from  <dig>  to  <dig>  ,  <dig>  to  <dig>   and  <dig>  to  <dig>   when comparing the 0% to the 50% missing value performance.
for each considered amount of missing data, d-pbtn outperformed depn ; results are comparable for aucpr, see table
<dig>
effect of noise
we next tested the effect of noise on inference performance, see figure
1c for aucroc, additional file
3: figure s <dig> for aucpr, and the summary in table
<dig>  to simulate noise, we randomly changed individual measurements in sn <dig> from  <dig> to  <dig> and vice versa.
network inference was performed on the resulting data , with parameters as described above, and resulting in comparable running times.
all four approaches showed decreasing performance for increasing levels of noise, with median aucroc decreasing from  <dig>   to  <dig>   for d-pbtn, from  <dig>  to  <dig>  for depn, from  <dig>  to  <dig>  for pbtn, and from  <dig>  to  <dig>  for bdagl.
d-pbtn was the best performing method using aucroc as performance measure, whereas pbtn was not better than guessing, with median aucroc =  <dig>  already at 16% noise.
in contrast to this, we observed good aucpr performance of pbtn, with aucpr =  <dig>  at all noise levels.
directly comparing d-pbtn with the non-time-series version pbtn, we observed that aucroc results for pbtn were inferior to d-pbtn at all noise levels, whereas pbtn showed better aucpr on data with higher levels of noise.
this is a very intersting observation, which we discuss further in the discussion section below.
bdagl showed very robust performance for increasing noise levels, consistently yielding low, but stable aucroc ≈  <dig> , but with very low aucpr ≈  <dig> .
we observed a further performance breakdown for bdagl when noise levels increased beyond 50% .
finally, comparing d-pbtn with depn, even with 50% of the data wrong, both approaches are significantly better than random guessing , and d-pbtn outperforms depn for all noise levels ; results for aucpr were qualitatively equivalent, see additional file
3: figure s <dig> and table
<dig>
overexpression experiments
as a further test of the method, we assessed the use of overexpression data with d-pbtn, using the sn1a data set.
model parameters were chosen as for the knockdown experiments above.
results were comparable to results obtained using the knockdown data, with aucroc =  <dig>   and a aucpr =  <dig>   for d-pbtn.
for comparison, we also used depn with the overexpression data.
this resulted in an aucroc =  <dig>  and aucpr =  <dig> , both inferior to the d-pbtn results.
we note that neither bdagl nor pbtn in its original implementation are able to handle overexpression data, and could therefore not be included in the comparison.
effect of network topology and model assumptions
data sn <dig> were simulated using the model underlying our approach.
results may thus favor d-pbtn and pbtn.
furthermore, artificial topologies may not be representative of real biological networks.
we therefore sampled "real" subgraphs from the kegg database, and simulated data using a different model .
the obtained data set sn <dig> was used to assess performance of d-pbtn, pbtn, bdagl and depn.
d-pbtn prior and mcmc and pbtn parameters were chosen as above, except for the parameter p of the binomial prior on τ, which was set to p =  <dig> , resulting in a slightly smaller expected value for τ. results for bdagl and depn were obtained as described above.
figure
1d and additional file
4: figure s <dig> show aucroc and aucpr results obtained, respectively, and all results are summarized in table
<dig>  median aucroc values achieved were  <dig>  ,  <dig>  ,  <dig>   and  <dig>  ; median aucpr values were  <dig>  ,  <dig>  ,  <dig>   and  <dig>  .
aucroc results for d-pbtn, pbtn and depn are statistically significantly better than guessing .
d-pbtn outperformed depn both with respect to aucroc and aucpr .
we observed that depn had particular difficulties when measured time points were far apart in time, with several unobserved intermediate steps between them.
d-pbtn, which in contrast to depn were specifically designed for time course data, did not exhibit this problem.
notably, the non time-series method pbtn showed only marginally inferior aucroc than d-pbtn , but superior performance in terms of median aucpr .
however, in contrast to d-pbtn, pbtn results are characterized by high variability
, iqr
), and mean aucpr values are very similar between the two approaches .
a statistical test comparing the aucpr values from the  <dig> replicate runs of the two approaches shows no significant difference in performance measured by aucpr .
we next evaluated the effect of negative feedback loops on network inference, using data set sn <dig>  this is a difficult problem, as the feedback can only be inferred from the temporal evolution of the network.
parameters for markov chains and prior were set as above, using  <dig> subpopulations with  <dig> chains each in sampling.
this resulted in a running time of  <dig> hours .
using the full data set, d-pbtn achieved a median aucroc of  <dig>  and median aucpr of  <dig> ; depns achieved a median aucroc of  <dig>  and median aucpr of  <dig> .
median aucroc performance values of pbtn and bdagl were both  <dig> , and were superior to depn, but inferior to d-pbtn.
in line with results on data set sn <dig>  we observed good aucpr performance of the non time-series method pbtn , however, in contrast to d-pbtn, pbtn could not infer the feedback loop between nodes  <dig> and  <dig>
additional file
5: figure s4a shows the distribution of edge weights obtained with d-pbtn, indicating that the posterior distribution  is unimodal.
using only observations of proteins  <dig> and  <dig>  unique networks can no longer be recovered .
probabilities for individual network topologies or even edges can be computed based on the sampled points, and can then be used e.g.
to design experiments to resolve the true underlying network .figure  <dig>
inference of negative feedback cycles.
the figure shows reconstructed network topologies on simulated data sn <dig>  using the full data set  and partial data with observations only for proteins  <dig> and  <dig> .
in the latter case, no unique network can be identified anymore.
clustering of sampled w groups topologies according to structural similarity.
only one cluster arises for the complete data , whereas four major different clusters arise on partial data .
shown are median networks from the four clusters, together with cluster probabilities computed as the fraction of networks in the cluster to the total number of sampled networks.
edge thickness indicates strength of support for the edge within the cluster , black lines indicate activations, red lines inhibitions.
robustness with respect to model hyperparameters
since our model contains parameters to be set by the user, we next analyzed how robustly networks were inferred for varying model hyperparameters, using data set sn <dig>  we varied γ in equation  by up to ±50%.
for changes of γ by the maximum change of ±50% tested, the change observed in inference performance was a decrease of  <dig> % for aucroc, and  <dig> % for aucpr .
of note, even for a change by 50%, d-pbtn still performed significantly better than random guessing .
we furthermore studied the influence of changing the prior hyperparameters q and s of the prior p ), changing them simultaneously by up to ±20% , or individually by up to 50% .
this resulted in a maximum change of  <dig> % in aucroc and  <dig> % in aucpr for individual changes of q by 50%, decreasing aucroc to  <dig>  and aucpr to  <dig> , and changes by  <dig> % and  <dig> % for aucroc and aucpr, respectively, for changes of s by 50%.
we then tested the influence of changing the parameters n and p of the binomial prior on τ by 50%, observing only a relatively minor impact of these parameters on performance .
a more significant influence was observed for the parameters r <dig> and s <dig> of the negative gamma prior on w <dig>  which influences the baseline activity level of unregulated nodes, see additional file
8: figure s <dig>  overall, robustness analysis indicates that inference is reasonably stable for varying hyperparameters; still, some care should be exercised in choosing values.
inferring erbb-mediated signal transduction
to test d-pbtn on real data, we used publicly available data from fröhlich et al.
<cit> , regarding  <dig> proteins involved in erbb signaling.
the data comprises  <dig> replicate measurements for  <dig> of the  <dig> proteins, measured after  <dig> different knock-downs using rna interference, including  <dig> combinatorial knockdowns.
measurements were taken at two different time points, directly before stimulation with egf, and after  <dig> hours of stimulation, using reverse phase protein arrays.
we discretized the normalized data using the midpoint of the negative and positive control medians as discretization threshold.
network inference was performed in  <dig> subpopulations with  <dig> chains each, sampling over  <dig> steps with a burn-in of  <dig>  steps.
parameters q and s were set to  <dig>  and  <dig>  respectively, corresponding to a fairly strict regularization.
we set the stochasticity parameter γ to γ =  <dig>  to obtain a slightly more deterministic behavior of the model than in case of the simulation study.
rate and shape parameters of the gamma prior p were set to r0 =  <dig> s0 =  <dig>  keeping unregulated nodes in the "off" state, and parameters of the binomial prior p were set to n =  <dig> and p =  <dig> , allowing for a relatively large number of model steps between the experimental measurements.
computation time took approximately  <dig> minutes, or roughly  <dig>  hours.
clustering of the  <dig> sampled weight vectors w shows only a single cluster with high probability, indicating a unimodal posterior.
to obtain a single network for visualization purposes and further biological discussion of inferred edges, we used a threshold δ =  <dig> × max|wi,j| for network discretization, including only edges with sample median
in the further evaluation - in words, edges were only included if the median of all points sampled for the corresponding edge was at least 65% of the maximum edge value sampled.
this results in a network with comparable number of edges as the gold standard network shown in additional file
9: figure s <dig>  constituting a tradeoff between false positive and false negative edges.
we compared results obtained using d-pbtn with the non-time series version pbtn
<cit> , as well as depn.
results for depn were taken from the original publication
<cit> , the authors do not report auc roc and pr values for their method.
figure
3a shows the resulting network using d-pbtn, and table
<dig> summarizes the results obtained using d-pbtn, pbtn and depn.
notably, without inclusion of any additional prior knowledge, both d-pbtn and depn show similar specificity and accuracy, but d-pbtn is both more sensitive and has higher precision than depn.
the dynamic version d-pbtn outperform pbtn in all performance measures used, showing the additional value of the temporal information in the data when using this method.figure  <dig>
reconstructed erbb-mediated signaling network.
the plots show inferred network topologies for reverse phase protein array data acquired after stimulation of erbb signaling and knockdowns of proteins involved in erbb signaling.
true positive edges, using high-confidence interactions in the string database as reference, are shown as solid black lines, newly predicted edges are shown as orange lines.
true and false negatives are not depicted for the sake of readability of the plot.
edge thickness indicates strength of support for the edge .
a: network inferred without using any prior knowledge, b: network inferred using the literature interactions reported by
<cit>  as prior knowledge for network inference.
performance comparison erbb signaling network
auc
roc
auc
pr
achieved sensitivity, specificity, accuracy and precision for the inference of erbb signaling
<cit> .
shown are results of d-pbtn, pbtn and depn, with and without prior information .
the last row shows results of the prior network alone.
network inference was assessed using string as gold standard.
results for depn were taken from froehlich et al.
<cit> , authors do not report auc
roc or auc
pr values.
the inferred network  is shown in figure
3a, and displays significant cross-talk in erbb-signaling.
at the surface receptor level, d-pbtn predicts cross-talk between egfr, esr and igfr, as recently confirmed experimentally
<cit> .
we furthermore predict a direct interaction between akt and mapk, linking the corresponding pathways.
in fact, involvement of mapk signaling in akt activation was recently reported
<cit> .
at the transcription factor level, d-pbtn predicts an interaction between cdk <dig> and myc, in line with observed expression changes of cdk <dig> and cdk <dig> following myc degradation
<cit> .
inference quality improves dramatically when prior knowledge is provided.
we used the literature network presented in fröhlich et al.
<cit>  as prior knowledge, with q =  <dig> and s =  <dig>  and μi,j = + <dig> in case of an activation, - <dig> in case of an inhibition, or  <dig> in case of no regulation in the prior network.
we note that this prior network alone yields superior results than d-pbtn without any prior knowledge , emphasising the importance of using a network prior on this data set.
network inference was performed as above, resulting in the network depicted in figure
3b.
obtained results are superior to both, using prior knowledge alone, as well as using only the experimental data .
to assess robustnes of these inference results with respect to model hyperparameters, we furthermore performed a sensitivity analysis by modifying model hyperparameters individually by up to ±50%, and evaluated changes in resulting auc roc and auc pr, compare figure
4a and b. this analysis indicated that results are relatively stable, with the most influential parameter being the stochasticity parameter γ.figure  <dig>
hyperparameter sensitivity analysis erbb network.
the figure shows sensitivity of the auc
roc  and auc
pr  of the erbb network  with respect to changes in model hyperparameters.
model parameters μ, s, q, γ, r
<dig> and s
<dig> where changed by up to ±50%, and inference was repeated for each parameter choice.
shown are resulting values for the area under the receiver operator characteristic  and the area under the precision-recall  curves.
the analysis indicates that inference results are reasonably robust to changes in model hyperparameters, with the stochasticity parameter γ being the most critical parameter for inference performance.
we evaluated predictions using ipa , and could find additional evidence for most of the predicted edges: d-pbtn predicts regulations of igfr, mapk and cdk <dig> by erbb <dig>  indeed, binding of erbb <dig> and igfr has been shown in lysate from skbr <dig> cells resistant to trastuzumab by immunoprecipitation
<cit> .
similarly, there is experimental evidence for an activation of mapk by erbb <dig>  as we predicted
<cit> .
the inferred activation of cdk <dig> by erbb <dig> to our knowledge has not been reported before, and may warrant further experiments to confirm or falsify this prediction.
some evidence also exists for the predicted strong activation of myc by esr.
both proteins have been shown by affinity chromatography to bind to one another
<cit> , and blocking of esr <dig> was recently shown to decrease myc expression
<cit> .
similarly, activation of akt increases ccnd and myc expression
<cit> , and a regulation of ccne by myc
<cit>  and the activation of ccnd by mapk
<cit>  have been demonstrated experimentally.
interestingly, d-pbtn predicts an activation of myc and ccne by igfr.
both interactions to our knowledge have not been previously reported, and may be interesting candidates for experimental validation.
discussion
in this manuscript, we present d-pbtn, a novel method to infer signal-transduction networks from time course perturbation data.
we evaluate d-pbtn in an extensive simulation study, and demonstrate its application on experimental data regarding erbb signaling.
the simulation study shows stable inference results in light of missing data and noise, and reasonable robustness with respect to model hyperparameters.
on the erbb application, we have run the inference with and without a network prior, and provide both results for comparison in figure
<dig>  while the overall topology is similar whether or not a prior network is used, there are significant differences in individual edges.
the quality of network inference with a network prior clearly depends on the quality of the prior used as well as on the amount and quality of experimental data available, we therefore recommend to carefully evaluate the effect a given network prior has.
in our example, there is experimental evidence for many of the predicted interactions, both from the run with and without network prior, directly confirming the quality of the inference.
furthermore, in addition to confirming known interactions, d-pbtn predicts several novel edges, which may be interesting candidates for experimental validation.
running time remains a major concern with d-pbtn, mainly due to the time required for sampling from the posterior distribution.
each single sampling step requires running time o, where k is the number of perturbations, t is the number of  time steps, and n is the number of nodes in the network.
this can be easily seen from the likelihood function , which essentially iterates over all time points, perturbations and nodes in the network, and computes equation  in each step, which again iterates over all nodes.
running time for the likelihood dominates the time required for the prior distribution, hence the total running time for the posterior is o. unfortunately, the shape and complexity of the likelihood function, and thus the number of sampling steps required to sample from the posterior p , also depends on k, t and n in a nontrivial way, and the choice of required sampling steps  is a nontrivial problem.
clearly running time increases linearly with the number of sampling steps, but how exactly the minimum number of sampling steps required depends on k, t and n is unclear.
while the present implementation of d-pbtn is a non-parallel implementation, evolutionary mcmc is straightforward to parallelize, allowing for a very efficient parallel implementation of the sampler.
in practical terms, running time on the sn <dig> example  was approximately  <dig>  hours with the present, non-parallelized version of the algorithm, and computations much beyond the size of the erbb network with its  <dig> nodes appear computationally prohibitive without further parallelization.
an efficient parallel implementation on a larger compute cluster in principle should make d-pbtn suitable also for "larger" network inference problems with several dozen, possibly up to a few hundred nodes, but this is clearly not an approach that is suitable for networks with thousands of proteins, let alone proteome-wide models.
as with all bayesian approaches, prior distributions and prior hyperparameters impact results, and need to be carefully chosen.
reasonable values for some of the parameters can be estimated from biological expectations, for example the expected level of "connectedness" in a network can be used to estimate parameters q and s of the prior on w, or information about the speed of signal transduction in a network could be used to estimate hyperparameters for the prior on τ. however, still a significant amount of "gut feeling" and experience are required to get good estimates.
while our simulation study indicates reasonable stability with respect to parameter changes, a possible extension of our work is to employ empirical bayes’ approaches to set model hyperparameters
<cit> , or, granted enough experimental data and compute time are available, crossvalidation-schemes could be used for hyperparameter estimation.
the comparison between d-pbtn and the non time-series version pbtn on the different simulated data sets gave some very interesting results.
while d-pbtn outperformed pbtn on all data sets when aucroc was used as performance measure, pbtn yielded higher median aucpr than d-pbtn on datasets with high levels of noise, or data simulated using a different underlying mathematical model.
aucpr in contrast to aucroc is less influenced by an unequal distribution of class labels , the superior results of pbtn over d-pbtn under certain conditions are therefore of relevance for the choice of inference method.
it can be shown that if a roc curve for method a dominates a roc curve for method b, then also the corresponding pr curve for method a dominates the pr curve for method b, and vice versa
<cit> .
the inverted ranking of d-pbtn and pbtn in aucroc and aucpr implies that the roc curves  of d-pbtn and pbtn do not dominate one another.
rather, there are certain domains in the roc and pr plots where one or the other method becomes better, and correspondingly the curves for d-pbtn and pbtn intersect.
this in particular occured in our simulation study for the data with high levels of noise.
it is tempting to speculate that under conditions with many "wrong" data points, focusing on a stable steady state  may have advantages over trying to capture the full dynamics of the system, as d-pbtn does.
this point will require further research in the future, to determine what exactly determines which of the two approaches is superior under given conditions, and to provide guidelines for method selection.
on the other hand, if the interest is in feedback loops and cycles in the networks, pbtn should not be used, as this approach cannot infer cyclic networks.
interestingly, on the erbb application example, d-pbtn clearly outperformed pbtn both with respect to aucroc and aucpr.
CONCLUSIONS
in this manuscript, we have developed a new method to infer signal transduction networks from time course perturbation data.
based on relatively few time points after a large number of different perturbations, our method is able to reconstruct the underlying signaling network with high accuracy.
the mathematical approach we employ is based on dynamical bayesian networks, and assumes discrete states and time steps, at which signaling molecules are either "on" or "off".
the approach is therefore suitable to experiments where protein activation or protein expression is completely or almost completely switched off, such as with rnai knockdowns.
several related approaches exist.
our work extends pbtn by explicitly taking time into account.
in contrast to pbtn, the method can thereby infer feedback and feedforward loops.
provided sufficient time-resolved measurements after  knockdowns are available, a unique network can be inferred using d-pbtn.
this is in contrast to methods such as pbtn, depn and nested effects models , which only return equivalence classes or single networks out of an equivalence class.
in this regard, our method is comparable to dynamic depns and dynamic nems, which also explicitly take time course data into account.
however, d-nems have been developed for "effect" observations after knockdowns, which only provide very indirect information about the underlying network.
d-depns, on the other hand, require long time series and work most efficiently if a relatively small number of perturbations are carried out.
the typical experimental scenario encountered in practice is different, comprising short time courses after a large number of perturbations.
d-pbtn have explicitly been developed for this situation.
due to the bayesian approach pursued, prior biological knowledge can easily be integrated into d-pbtn inference, and marginalization over unobserved proteins or time points is straightforward and very efficient due to the likelihood simulation.
we employ an mcmc approach to sample from the posterior distribution of models given the experimental data, which allows it to compute probability distributions over alternative network topologies.
this is particularly useful if multiple models explain the data equally well, and can be used to design additional experiments to then distinguish further between high-scoring topologies.
using halite clustering, similar models can be grouped together and probabilities estimated from the mcmc samples.
overall, d-pbtn is a powerful approach for network inference when time-resolved measurements of the response of a dynamical system after perturbation are available.
under these conditions, the approach outperforms other state-of-the-art methods.
most approaches available in the literature consider either changes in steady state levels after perturbation, or the temporal dynamics of an unperturbed system for network inference.
d-pbtn makes use of both the dynamical aspects and the perturbation effects in network reconstruction.
we present an extensive simulation study to evaluate the approach, showing stable performance under varying conditions of noise and missing values in the data, and we demonstrate the application of d-pbtn to infer signaling networks in the erbb system.
here, our approach both reconfirms known interactions, as well as suggesting some novel edges as candidates for further experimental study.
availability and requirements
project name: dynamic probabilistic threshold networks
project home page:http://www.kaderali.org/software/dpbtn.html.
operating system: tested on microsoft windows
programming language: c#
other requirements: none
license: gnu gpl
any restrictions to use by non-academics: none
electronic supplementary material
additional file 1:
supplementary material.
additional file one contains supplementary text with additional information.
additional file 2:
supplementary figure s <dig>  effect of the amount of unobserved proteins on inference performance .
shown is the distribution of the area under the precision recall curve  for  <dig> replicates of simulated data sets, over the fraction of unobserved proteins, for d-pbtn, depn, bdagl and pbtn on dataset sn <dig>  the dashed line shows expected auc
pr results for random guessing.
additional file 3:
supplementary figure s <dig>  effect of noise on inference performance .
the figure shows the distribution of auc
pr values obtained for different levels of noise on sn <dig>  for d-pbtn, depn, bdagl and pbtn.
the dashed line shows expected auc
pr results for random guessing.
additional file 4:
supplementary figure s <dig>  inference performance on the reconstruction the kegg networks.
this figure shows performance of the d-pbtn, depn bdagl and pbtnapproaches on subnetworks taken from the kegg database, with data simulated as described in methods.
shown are auc
pr values of  <dig> simulated data sets, generated from ten different kegg subnetworks.
the dashed line shows expected auc
pr results for random guessing.
additional file 5:
supplementary figure s <dig>  inference of negative feedback cycles: bean-plot showing the distributions of inferred edge weights.
a unimodal distribution is observed with the full data , whereas no unique network can be identified when only observations of two downstream proteins are given.
additional file 6:
supplementary figure s <dig>  effect of model parameters on performance: the figure shows the effect of varying the stochasticity parameter on the auc values .
network inference has been done using data simulated from simulated network  <dig>  auc
roc values are shown as red squares, auc
pr values are shown as blue circles.
additional file 7:
supplementary figure s <dig>  effect of model parameters on performance: the figure shows the effect of varying the rate and shape parameter q and s of the prior p on the auc values.
network inference has been done using data simulated from simulated network  <dig>  panels a and b: effect of simultaneous changes of q and s by up to 20% on auc
roc and auc
pr.
panels b and c: effect of changing parameters q and s individually by up to 50%.
additional file 8:
supplementary figure s <dig>  effect of model parameters on performance: the figure shows the effect of varying the rate and shape parameters r
<dig>  and s
<dig>  of the prior on the bias weight w
<dig> by up to 50%.
network inference has been done using data simulated from simulated network  <dig>  auc
roc values are shown as red squares, auc
pr values are shown as blue circles.
additional file 9:
supplementary figure s <dig>  experimental reference network extracted from string for erbb network inference.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
lk and nak designed the research and methodology.
nak and lk wrote the software and analyzed the data.
lk and nak wrote the paper.
both authors read and approved of the final version of the manuscript.
we acknowledge funding from the german ministry of education and research  via systec/fancy  and gerontosys/agenet , and from the european union 7th framework program .

 the availability of high throughput experimental platforms has transformed molecular biology into a data-rich science. however, the development of computational and mathematical tools to extract and interpret the wealth of information hidden in these data is lagging behind. for example, genome wide rna interference screens have enabled the phenotypic characterization of genes at an unprecedented scale in living cells. however, the interpretation of such data and the placement of hits in their functional and temporal context in cellular pathways remains a major challenge. machine learning can address many of the questions arising from the interpretation of large scale molecular biological data, and has become an important tool of bioinformatics and systems biology. automatic network inference deals with the reconstruction of regulatory or signaling networks directly from experimental data using statistical or machine learning approaches, a field that has received significant attention in the last decade. there is a wide variety of different approaches available that can be used to infer genetic regulatory or signal transduction networks from experimental data. approaches employed include bayesian networks and dynamic bayesian networks, boolean models, auto-regressive models, correlation-based and mutual-information based models, clustering techniques, differential equation models, and others. these methods differ in the level of detail at which they reconstruct networks, and in their underlying assumptions and data requirements. some produce an undirected graph, where edges do not indicate which gene or protein in a connected pair is the activator, and which is the activated gene or protein. others specify the regulator with directed edges, and a few label the edges with kinetic parameters. there are a number of successful applications of network inference approaches to elucidate cellular signaling pathways, including meta-approaches integrating different methods. among the successful applications, for example, sachs et al. used bayesian networks to reconstruct cellular protein signaling networks from protein phosphorylation measurements. nelander et al. used a model based on nonlinear differential equations to infer signaling networks in cancer from combinatorial drug perturbation data. eduati et al. demonstrate the integration of literature-knowledge into data driven approaches, and show a successful application to a signaling network related to growth-signaling and inflammation. ciaccio et al. used bayesian networks and two different mutual-information based approaches to infer signaling networks downstream of the egf receptor. hill et al. use dynamic bayesian networks to study signaling in a cancer cell line. other approaches include nested effects models, deterministic effects propagation networks or probabilistic boolean threshold networks, and have been applied, for example, to reconstruct signaling networks in the erbb pathway, or in the innate immune response to infection. network reconstruction can be performed from observational data alone. however, the quality of the reconstruction increases substantially if experimental perturbations followed by observations of the system’s dynamic response are available. surprisingly, while there are many approaches available for time course data, and several approaches for perturbation data, there are only few methods available that can handle both time course data and perturbations at the same time. exceptions are, for example, the differential equation approach presented by nelander et al., dynamic nested effects models or dynamic deterministic effects propagation models. among the stochastic approaches, dnems rely on high-dimensional, indirect readouts of rather qualitative knockdown "effects", such as microarrays performed at different time points after every knockdown, or multidimensional features derived from live cell imaging. such data are often not available, and also provide only indirect information about the underlying signaling pathway. fröhlich et al. proposed deterministic effects propagation networks to reconstruct networks from perturbation data with direct observation of involved proteins, measured e.g. using reverse phase protein arrays. however, depns treat each time point as an independent measurement and do not model the time dependent behavior of the system explicitly. bender et al. proposed dynamic depns, to take the availability of longitudinal data as well as inhibitory interactions in a network into account. d-depns use the viterbi algorithm to identify the state transitions of a hidden markov model from the time course data, with states corresponding to combinations of activities of nodes in the network. a likelihood function is defined to score network models given the estimated state transitions, and the network space is then searched using a genetic algorithm. due to this underlying procedure, d-depns require relatively long time series, and have substantial running time requirements if large numbers of different perturbations involving most or all of the proteins in a network are performed. we here focus on the problem of reconstructing signaling networks from short time series with only two or three time points, after a large number of different, possibly combinatorial perturbations, targeting most or all of the genes in the network. such data arise, for example, if rnai experiments are combined with protein array measurements at few time points after the perturbation, or in time-resolved mass spectrometric assays under different conditions. we propose a new method for the identification of kinetic models of signaling networks from such data, dynamic probabilistic boolean threshold networks, which can treat multiple, combinatorial interventions as well as incomplete observations. through the use of a fairly simple, discrete-state model where proteins are either active or inactive, our method is applicable to qualitative data and does not require detailed quantitative measurements. furthermore, through the bayesian framework employed, we can easily handle noisy or missing data, and using mcmc sampling, can analyze full posterior distributions over model parameters, yielding information about different, alternative network topologies that are consistent with the experimental data. our work builds upon probabilistic boolean threshold networks regarding the underlying dynamic network model, but extends the method in the following aspects: we fully exploit the information from time course data, by explicitely integrating time into the model likelihood. in contrast, pbtns in their original implementation only exploit observational data taken at a single time point, usually at steady state. pbtns are limited to relatively small networks of around nodes. through a more efficient sampling method, the approach presented here can be used for larger networks, effectively more than doubling the feasible network size. pbtns were implemented with a focus on data with only single downstream observations at a single time point, and therefore similarly to other bayesian network based approaches, do not allow feedback loops in the network. the approach presented here can handle feedback loops in the network. d-pbtns can handle overexpression and knockdown data, in contrast to pbtns, which are limited to knockdown data. our approach uses discrete states of the proteins in the network, which can be either "on" or "off". in terms of the biological interpretation, this could be phosphorylation networks, in which case edges would correspond to phosphorylation or dephosphorylation events. alternatively, assuming that proteins are "present" or "absent", edges in the model could also be interpreted as transcriptional regulation. correspondingly, nodes can represent "activated" protein levels or total protein abundance depending on context. we evaluate the performance of our method on simulated data, and study its behavior on networks of different sizes, with different amounts of available data and different levels of noise, and compare results against pbtn, depn and bdagl, a dynamic bayesian network based approach. dynamic depns could not be used for the comparison, as they were unable to handle the short time courses and large number of perturbations. we show inference results for erbb signaling in breast cancer using depn, pbtn and d-pbtn, demonstrating superior performance of d-pbtn, and showing cross-talk between different branches of the erbb pathway. methods mathematical model we model a signaling network by a weighted, directed graph. cycles are permitted. the node set corresponds to proteins; an edge from node vi to vj represents a regulation of vj by vi, such as activation or deactivation via phosphorylation. the strength and type of an interaction is specified by, with wi,j >  for activations, wi,j <  for inactivations and wi,j:=  if. given n, the number of nodes in the network, we can thus describe the network topology by its n × n weighted adjacency matrix. we furthermore consider each node vi a boolean variable. model time is assumed discrete, and all nodes are updated simultaneously between two time steps. the state of node vi at model time t +  is described by a stochastic function of the states of all nodes at model time t: the relationship between a target protein and its regulators is hence modeled by a sigmoid function, an ansatz for nonlinear systems with saturation phenomena. equation activates or deactivates a target protein depending on the weighted sum of incoming regulations. the level of stochasticity in this process is defined by the positive parameter γ, whereas w i defines a basal activation probability of vi in the absence of any regulatory effects. in case of interventions, the state of a node v targeted by perturbation is no longer subject to the stochastic dynamics of the system, instead, its value is fixed through the intervention. we write, k ∈ { …,k}, for the set of nodes targeted simultaneously in perturbation experiment k, and assume that a knockdown of amounts to fixing all nodes to the "off" state at all model time points after the knockdown. similarly, overexpression experiments can be simulated by fixing the affected nodes vi to "on". let us now assume that we have observed experimentally the states of the nodes vi at t different time points,. we note in passing that this normally requires a discretization step, as experimental data are typically measured continuously – see the section on the erbb data set below. let us for now assume that the set of experimental observations consists of measurements of n discrete protein states at t time points after k different perturbations. we write for the observed activity of node vi at the real time point after knockdown of node set. to map real time to model time t, we introduce a parameter vector to denote the number of model time steps that correspond to the real time elapsed between the different experimental time points, that is, denotes the number of model time steps according to equation required to transition from real time point to time point. assuming a one-to-one correspondence between model time and real time, that is, assuming that for all, we can define a likelihood by calculating the product of according to equation over all nodes v, time points t and perturbations k: if τt >  one or more intermediate model time steps must be made to transition from experimental time point to experimental time point. as no experimental observations are available for the intermediate time steps, we must marginalize over unobserved time steps in the evaluation of the likelihood. writing for the full sequence of model steps, where the sequence of intermediate steps are unobserved, equation becomes where comprises the product over all individual nodes as in equation, and where the sum is over all possible combinations of values for the intermediate model steps. missing values in the experimental data can be treated similarly – the likelihood can then be computed by integrating equation over all possible combinations of values for unobserved. we note that pbtns are based on a closely related likelihood function, and employ the same underlying model for the relation between nodes, given in equation. however, pbtns differ in a key aspect from the framework presented here: pbtns do not allow multiple time points in their likelihood, but rather assume that measurements are made at a single time point after knockdown, usually taken at steady state. inclusion of prior knowledge in many biological settings, available data are insufficient to unambiguously reconstruct the underlying network. in such situations, strict regularization of the objective function using for example maximum parsimony, or the inclusion of additional prior biological knowledge can help. both can be done via bayes’ theorem using prior distributions on the model parameters Ω. the model posterior is where is the likelihood function as defined above, p is an adequately defined prior distribution, and is a constant normalizing factor that does not depend on Ω and can be neglected in maximization or sampling. the prior can be written as, assuming that the different types of model parameters are independent. then, describes our belief about the correct topology, prior to seeing any data, and p describes our belief about the speed of signal transduction in the network. assuming independence of parameters, the full prior can be written as: in the network inference setting, it is unlikely that the true underlying biological network is densely connected. we rather expect a sparse network, where most of the nodes have only few other nodes acting on them. this can be mathematically encoded in a prior on wi,j with positive shape and rate parameters q and s. note that, after taking the negative logarithm and dropping the normalizing factor 1/, this is equivalent to regularization using the lq-norm, and corresponds to a laplacian prior if q =  or a normal prior if q =  and the parameter s determines the width of the distribution. the parameters μi,j can be used to encode additional knowledge on specific edges, by setting individual μi,j to nonzero values. if no prior knowledge is available, μi,j defaults to we note that the μi,j can be chosen continuously depending on the expected strength of the effect and certainty in its presence. we expect signaling molecules to be "off" in the absence of any signal, and therefore use independent negative gamma priors for p, with rate and shape parameters r and s respectively. we hence assume w i to be non-positive. this assumption can easily be replaced if signaling molecules should also be allowed to be "on" in the absence of signal, for example by using a gamma prior on the absolute value |w i|, or a normal prior on w i. lastly, as model time is discrete and necessarily positive, we use a binomial prior p on τi, with parameters n and p. network inference we could now maximize to find the most probable network Ω underlying the data. this is reasonable if the posterior is unimodal and has a sharp, narrow optimum. however, since the optimization problem encountered for real data is typically underdetermined, several alternative models Ω often explain the data equally well, and marginal distributions for parameters show wide peaks and large confidence intervals. we therefore use markov chain monte carlo sampling to fully characterize the posterior. the main difficulty in doing this lies in the calculation of the likelihood, which due to the required marginalization over unobserved nodes and time points quickly becomes very involved and time-consuming, compare equation. we therefore approximate by simulating the dynamics of the model with parameters Ω for each knockdown, and then compare the simulation results mk at each time point t with the experimental data dk. if this is repeated a large number of times, we can approximate by the number of times we get back in the simulation runs. by combining this with mcmc sampling, a very efficient approach for evaluation of the posterior distribution arises, that does not require an explicit computation of the likelihood. in contrast to pbtn, which employed a metropolis-hastings based sampler, we have now implemented a far more efficient sampler in d-pbtn. we employ distributed evolutionary monte carlo, a sampler that combines features of a distributed genetic algorithm with mcmc sampling. in brief, demc starts with a population of m markov chains, which are grouped into g subpopulations. each individual chain describes a specific network Ω. updates within each subpopulation are done using the genetic operators mutation and cross-over, and migration allows individual solutions to move between subpopulations. thereafter, each chain is scored by initializing each node in the network with the values experimentally observed at time and then simulating using equation. this is followed by a metropolis-hastings step to accept or reject the new state, to ensure ergodicity of the chains. a detailed description of the procedure is given in additional file data simulation since for real experimental data typically no "gold standard" network is available to assess results, we used simulated data to evaluate d-pbtn. this allows it to systematically alter properties of the data, such as the inherent level of noise or amount of missing values, and to directly assess the performance impact this has. we simulated data in three different ways: simulated network is a 7-node feedforward network, shown in figure 1a. weights wi,j were set to for edges shown in the figure, and otherwise. baseline weights were set to w i = - to have a high probability for an unregulated protein to be in the off state. data simulation for this network was performed assuming that all proteins are in the off state initially, except for the receptor stochastic signaling is then simulated using equation, updating all nodes simultaneously, and using γ =  knockdowns were simulated by fixing targeted proteins in the off state. we simulated knockdowns of all individual proteins, and two combinatorial knockdowns of proteins and. two time points were used for network inference, the first one immediately after knockdown and activation of the receptor; the second time point after simulation steps, allowing for sufficient time for the signal to propagate through the whole network.figure evaluation on simulated data. a: the panel shows network topologies used to simulate data. simulated network is a simple feedforward network, whereas simulated network contains a negative feedback loop. b: network reconstruction was performed for sn changing the number of unobserved proteins. shown is the distribution of the area under the roc curve of replicates of simulated data sets, over the fraction of unobserved proteins. from left to right: d-pbtn, depn, bdagl, pbtn. the dashed line at auc roc =  shows expected results for random guessing. c: the panel shows the distribution of auc roc values obtained for different levels of noise on sn noise is introduced by switching the state of the indicated fraction of proteins in the "measured" data, thus introducing errors in the data. d: this panel shows performance of d-pbtn, depn, bdagl and pbtn on networks sampled from kegg, with data simulated as described in methods. shown are auc roc values of simulated data sets, generated from ten different kegg subnetworks. using the same network topology and parameters, we furthermore simulated overexpression experiments. overexpressed nodes were fixed to the on state, while we assumed all other nodes to be inactive initially. we simulated overexpression experiments of all individual nodes, as well as combinatorial overexpression of nodes and nodes, again by running the simulation over time steps as above. simulated network extends sn by a feedback loop, see figure 1a. data for sn was simulated using the model proposed by fröhlich et al.. in contrast to sn signaling is deterministic, and noise arises only at the measurement stage. we initialized all nodes as above and assumed that a node deterministically becomes active if there are more activating than inhibiting proteins among its parents. measurements were simulated by sampling from a normal distribution with -distributed mean for active nodes, and -distributed mean if the node was inactive. the variances of the normal distributions were drawn from an inverse scaled chi-squared distribution, as suggested by fröhlich et al. based on their experimental observations. measurements were simulated immediately after the knockdown and activation, and after four and six simulation steps. simulated network comprises a set of ten different networks randomly extracted from the kegg database, each containing seven connected proteins, taking only protein-protein interactions into account. we simulated time-course data for time points as described for sn for each of the ten subnetworks. for each subnetwork, in addition to the starting point, two further randomly selected time points were included into the final data set used for network reconstruction. reported performance results are average performances over all ten subnetworks, avoiding biasing of results towards a single topology or time interval. implementation and performance evaluation we implemented d-pbtn in c♯, and evaluated the approach on both simulated and real experimental data. source code and additional information is available from our website at http://www.kaderali.org/software.html. convergence of the markov chains was assessed using the methods proposed by brooks and gelman, results are given in additional file sampled networks were aggregated using halite clustering. halite is a density-based clustering method, that we use to identify clusters of similar networks from the markov chains. cluster representative networks were computed using the within-cluster median of wi,j. to quantify network reconstruction performance, we performed receiver operator characteristic and precision recall analysis. in brief, a threshold δ is used for discretization of edge values, where the median of sampled values for a particular edge wi,j is compared with δ. if, the edge ei,j is assumed present in the network,, and absent otherwise,. the resulting network with edge set is then compared against the gold standard network, and sensitivity, specificity and precision are computed for given δ. this is then repeated by varying δ, and sensitivity is plotted over specificity for different δ in a roc plot, and precision over sensitivity in a pr plot. finally, both roc and pr curves can be summarized by computing the area under the curve as single numbers summarizing performance over a wide range of threshold values δ. the advantage of this approach is that it is not necessary to pick a specific threshold δ for the discretization, which may be difficult to do as this depends on user preferences regarding the tradeoff between false positive and false negative edges, but the analysis summarizes performance over the full range of feasible thresholds δ. we first assessed performance of d-pbtn on simulated data, varying the amount of noise in the data, using different amounts of simulated data and different models for data simulation, and evaluated stability of results for changing model hyperparameters. we then assessed performance on publicly available data regarding signaling in the erbb network. inferred networks on real data were assessed using string as reference, using only edges with at least 70% confidence in string. since interactions given in string are undirected and unsigned, we disregarded directional information and edge sign. we compare performance of our approach with depn, the non-time-series version of our approach, the bayesian directed acyclic graph learning tool developed by eaton and murphy, and with random guessing. notably, d-depn could not be used on the simulated data due to the nature of the short time courses used. results for pbtn were obtained using a modified version of the original pbtn implementation that optimizes the posterior distribution instead of sampling, thus allowing for a more efficient computation on networks with a unimodal posterior. this compromise is necessary, as the original sampling-version of pbtn is too slow to allow a rigorous evaluation on networks of the size used here. time points were used as independent measurements in pbtn. to obtain results for the depns, we used greedy hill-climbing and bootstrapping with bootstrap samples, as proposed by the authors of the method. the bootstrap samples were used to obtain weights on the edges, which were subjected to roc and pr analysis. we note that depns operate on equivalence classes of graphs, and can not distinguish between graphs within an equivalence class. we therefore used the transitively closed network as a representative of the full equivalence class to evaluate the depn results. this is in contrast to the other methods, which in principle are able to determine a unique network, provided sufficient data is given, and which we compared directly to the gold standard network. lastly, results for bdagl were obtained using the modified marginal likelihood method for perfect interventional data, with uniform prior.  as a first test of our method, we assessed the performance of d-pbtn on simulated data, without noise and with full observations of all nodes, using the network topology sn to compare performance of d-pbtn with other approaches, we used the non-time series version pbtn, depn and bdagl. for d-pbtn, shape and rate parameters of the prior distribution on w were set to q,s =  respectively, corresponding to a standard l regularization of the edge weights as the simplest conceivable setting. rate and shape parameters of the negative gamma prior on w i were set to r0 =  s0 =  resulting in a moderate deactivation of unregulated proteins, and a binomial prior with n =  p =  was used for τ, giving an expected value of τ of time steps. the stochasticity parameter γ was set to, corresponding to an average level of noise in the data. three sub-populations with four markov chains were used with steps each and a burn-in of steps. parameters for pbtn were set accordingly, using the two time points as independent replicate measurements. results for all four approaches are shown in table values shown are median values of the aucroc and aucpr measures over replicate simulations. with no noise and no missing values, d-pbtn shows superior performance over the other three approaches on the simulated data. notably, the comparison between d-pbtn and pbtn shows the added value of the temporal information, indicating that while a significant part of the information in the data comes from the knockdowns, even short time courses with just time points are of added value. bdagl shows weaker results than the other three approaches, both in terms of the aucroc and aucpr, possibly due to a significantly different underlying mathematical model.table performance comparison of d-pbtn, pbtn, depn and bdagl on simulated data auc roc auc pr network sn1 sn2 sn3 shown are achieved values of the area under the roc curve and the area under the pr curve. values shown are median values over iterations. inference was performed on the sn data, with data with different levels of noise and missing values, on the sn data including a negative feedback loop, and on the sn data set. note that bdagl cannot handle missing values. the upper part of the table shows the area under the roc curve, while the lower part of the table shows the area under the pr curve. effect of missing data we next assessed the effect of missing observations on network inference, using the sn data. results are shown in figure 1b for aucroc and additional file 2: figure s for aucpr, and are summarized in table note that bdagl could not be used in this comparison, as the method cannot handle missing values. we randomly selected 16%, 33% and 50% of the nodes, and completely removed the corresponding simulated observations for these nodes before carrying out the network inference. model parameters were set as described above. computation took approximately hours for d-pbtn on a ghz intel 64-bit processor. for comparison, runtime for depn was only around minutes. regarding quality of the inferred networks, d-pbtn shows superior performance over depn, both with respect to aucroc and aucpr. pbtn did not perform better than random guessing on all runs with missing values. both d-pbtn and depn show a decrease in the quality of the inferred network with increasing amounts of missing observations. in spite of this, even with 50% of missing data, both approaches still perform significantly better than guessing. median aucroc decreased from to, to and to when comparing the 0% to the 50% missing value performance. for each considered amount of missing data, d-pbtn outperformed depn; results are comparable for aucpr, see table effect of noise we next tested the effect of noise on inference performance, see figure 1c for aucroc, additional file 3: figure s for aucpr, and the summary in table to simulate noise, we randomly changed individual measurements in sn from to and vice versa. network inference was performed on the resulting data, with parameters as described above, and resulting in comparable running times. all four approaches showed decreasing performance for increasing levels of noise, with median aucroc decreasing from to for d-pbtn, from to for depn, from to for pbtn, and from to for bdagl. d-pbtn was the best performing method using aucroc as performance measure, whereas pbtn was not better than guessing, with median aucroc =  already at 16% noise. in contrast to this, we observed good aucpr performance of pbtn, with aucpr =  at all noise levels. directly comparing d-pbtn with the non-time-series version pbtn, we observed that aucroc results for pbtn were inferior to d-pbtn at all noise levels, whereas pbtn showed better aucpr on data with higher levels of noise. this is a very intersting observation, which we discuss further in the discussion section below. bdagl showed very robust performance for increasing noise levels, consistently yielding low, but stable aucroc ≈ , but with very low aucpr ≈ . we observed a further performance breakdown for bdagl when noise levels increased beyond 50%. finally, comparing d-pbtn with depn, even with 50% of the data wrong, both approaches are significantly better than random guessing, and d-pbtn outperforms depn for all noise levels; results for aucpr were qualitatively equivalent, see additional file 3: figure s and table overexpression experiments as a further test of the method, we assessed the use of overexpression data with d-pbtn, using the sn1a data set. model parameters were chosen as for the knockdown experiments above. results were comparable to results obtained using the knockdown data, with aucroc =  and a aucpr =  for d-pbtn. for comparison, we also used depn with the overexpression data. this resulted in an aucroc =  and aucpr = , both inferior to the d-pbtn results. we note that neither bdagl nor pbtn in its original implementation are able to handle overexpression data, and could therefore not be included in the comparison. effect of network topology and model assumptions data sn were simulated using the model underlying our approach. results may thus favor d-pbtn and pbtn. furthermore, artificial topologies may not be representative of real biological networks. we therefore sampled "real" subgraphs from the kegg database, and simulated data using a different model. the obtained data set sn was used to assess performance of d-pbtn, pbtn, bdagl and depn. d-pbtn prior and mcmc and pbtn parameters were chosen as above, except for the parameter p of the binomial prior on τ, which was set to p = , resulting in a slightly smaller expected value for τ. results for bdagl and depn were obtained as described above. figure 1d and additional file 4: figure s show aucroc and aucpr results obtained, respectively, and all results are summarized in table median aucroc values achieved were,, and; median aucpr values were,, and. aucroc results for d-pbtn, pbtn and depn are statistically significantly better than guessing. d-pbtn outperformed depn both with respect to aucroc and aucpr. we observed that depn had particular difficulties when measured time points were far apart in time, with several unobserved intermediate steps between them. d-pbtn, which in contrast to depn were specifically designed for time course data, did not exhibit this problem. notably, the non time-series method pbtn showed only marginally inferior aucroc than d-pbtn, but superior performance in terms of median aucpr. however, in contrast to d-pbtn, pbtn results are characterized by high variability, iqr ), and mean aucpr values are very similar between the two approaches. a statistical test comparing the aucpr values from the replicate runs of the two approaches shows no significant difference in performance measured by aucpr. we next evaluated the effect of negative feedback loops on network inference, using data set sn this is a difficult problem, as the feedback can only be inferred from the temporal evolution of the network. parameters for markov chains and prior were set as above, using subpopulations with chains each in sampling. this resulted in a running time of hours. using the full data set, d-pbtn achieved a median aucroc of and median aucpr of; depns achieved a median aucroc of and median aucpr of. median aucroc performance values of pbtn and bdagl were both, and were superior to depn, but inferior to d-pbtn. in line with results on data set sn we observed good aucpr performance of the non time-series method pbtn, however, in contrast to d-pbtn, pbtn could not infer the feedback loop between nodes and additional file 5: figure s4a shows the distribution of edge weights obtained with d-pbtn, indicating that the posterior distribution is unimodal. using only observations of proteins and unique networks can no longer be recovered. probabilities for individual network topologies or even edges can be computed based on the sampled points, and can then be used e.g. to design experiments to resolve the true underlying network.figure inference of negative feedback cycles. the figure shows reconstructed network topologies on simulated data sn using the full data set and partial data with observations only for proteins and. in the latter case, no unique network can be identified anymore. clustering of sampled w groups topologies according to structural similarity. only one cluster arises for the complete data, whereas four major different clusters arise on partial data. shown are median networks from the four clusters, together with cluster probabilities computed as the fraction of networks in the cluster to the total number of sampled networks. edge thickness indicates strength of support for the edge within the cluster, black lines indicate activations, red lines inhibitions. robustness with respect to model hyperparameters since our model contains parameters to be set by the user, we next analyzed how robustly networks were inferred for varying model hyperparameters, using data set sn we varied γ in equation by up to ±50%. for changes of γ by the maximum change of ±50% tested, the change observed in inference performance was a decrease of % for aucroc, and % for aucpr. of note, even for a change by 50%, d-pbtn still performed significantly better than random guessing. we furthermore studied the influence of changing the prior hyperparameters q and s of the prior p ), changing them simultaneously by up to ±20%, or individually by up to 50%. this resulted in a maximum change of % in aucroc and % in aucpr for individual changes of q by 50%, decreasing aucroc to and aucpr to, and changes by % and % for aucroc and aucpr, respectively, for changes of s by 50%. we then tested the influence of changing the parameters n and p of the binomial prior on τ by 50%, observing only a relatively minor impact of these parameters on performance. a more significant influence was observed for the parameters r and s of the negative gamma prior on w which influences the baseline activity level of unregulated nodes, see additional file 8: figure s overall, robustness analysis indicates that inference is reasonably stable for varying hyperparameters; still, some care should be exercised in choosing values. inferring erbb-mediated signal transduction to test d-pbtn on real data, we used publicly available data from fröhlich et al., regarding proteins involved in erbb signaling. the data comprises replicate measurements for of the proteins, measured after different knock-downs using rna interference, including combinatorial knockdowns. measurements were taken at two different time points, directly before stimulation with egf, and after hours of stimulation, using reverse phase protein arrays. we discretized the normalized data using the midpoint of the negative and positive control medians as discretization threshold. network inference was performed in subpopulations with chains each, sampling over steps with a burn-in of steps. parameters q and s were set to and respectively, corresponding to a fairly strict regularization. we set the stochasticity parameter γ to γ =  to obtain a slightly more deterministic behavior of the model than in case of the simulation study. rate and shape parameters of the gamma prior p were set to r0 =  s0 =  keeping unregulated nodes in the "off" state, and parameters of the binomial prior p were set to n =  and p = , allowing for a relatively large number of model steps between the experimental measurements. computation time took approximately minutes, or roughly hours. clustering of the sampled weight vectors w shows only a single cluster with high probability, indicating a unimodal posterior. to obtain a single network for visualization purposes and further biological discussion of inferred edges, we used a threshold δ =  × max|wi,j| for network discretization, including only edges with sample median in the further evaluation - in words, edges were only included if the median of all points sampled for the corresponding edge was at least 65% of the maximum edge value sampled. this results in a network with comparable number of edges as the gold standard network shown in additional file 9: figure s constituting a tradeoff between false positive and false negative edges. we compared results obtained using d-pbtn with the non-time series version pbtn, as well as depn. results for depn were taken from the original publication, the authors do not report auc roc and pr values for their method. figure 3a shows the resulting network using d-pbtn, and table summarizes the results obtained using d-pbtn, pbtn and depn. notably, without inclusion of any additional prior knowledge, both d-pbtn and depn show similar specificity and accuracy, but d-pbtn is both more sensitive and has higher precision than depn. the dynamic version d-pbtn outperform pbtn in all performance measures used, showing the additional value of the temporal information in the data when using this method.figure reconstructed erbb-mediated signaling network. the plots show inferred network topologies for reverse phase protein array data acquired after stimulation of erbb signaling and knockdowns of proteins involved in erbb signaling. true positive edges, using high-confidence interactions in the string database as reference, are shown as solid black lines, newly predicted edges are shown as orange lines. true and false negatives are not depicted for the sake of readability of the plot. edge thickness indicates strength of support for the edge. a: network inferred without using any prior knowledge, b: network inferred using the literature interactions reported by as prior knowledge for network inference. performance comparison erbb signaling network auc roc auc pr achieved sensitivity, specificity, accuracy and precision for the inference of erbb signaling. shown are results of d-pbtn, pbtn and depn, with and without prior information. the last row shows results of the prior network alone. network inference was assessed using string as gold standard. results for depn were taken from froehlich et al., authors do not report auc roc or auc pr values. the inferred network is shown in figure 3a, and displays significant cross-talk in erbb-signaling. at the surface receptor level, d-pbtn predicts cross-talk between egfr, esr and igfr, as recently confirmed experimentally. we furthermore predict a direct interaction between akt and mapk, linking the corresponding pathways. in fact, involvement of mapk signaling in akt activation was recently reported. at the transcription factor level, d-pbtn predicts an interaction between cdk and myc, in line with observed expression changes of cdk and cdk following myc degradation. inference quality improves dramatically when prior knowledge is provided. we used the literature network presented in fröhlich et al. as prior knowledge, with q =  and s =  and μi,j = + in case of an activation, - in case of an inhibition, or in case of no regulation in the prior network. we note that this prior network alone yields superior results than d-pbtn without any prior knowledge, emphasising the importance of using a network prior on this data set. network inference was performed as above, resulting in the network depicted in figure 3b. obtained results are superior to both, using prior knowledge alone, as well as using only the experimental data. to assess robustnes of these inference results with respect to model hyperparameters, we furthermore performed a sensitivity analysis by modifying model hyperparameters individually by up to ±50%, and evaluated changes in resulting auc roc and auc pr, compare figure 4a and b. this analysis indicated that results are relatively stable, with the most influential parameter being the stochasticity parameter γ.figure hyperparameter sensitivity analysis erbb network. the figure shows sensitivity of the auc roc and auc pr of the erbb network with respect to changes in model hyperparameters. model parameters μ, s, q, γ, r and s where changed by up to ±50%, and inference was repeated for each parameter choice. shown are resulting values for the area under the receiver operator characteristic and the area under the precision-recall curves. the analysis indicates that inference results are reasonably robust to changes in model hyperparameters, with the stochasticity parameter γ being the most critical parameter for inference performance. we evaluated predictions using ipa, and could find additional evidence for most of the predicted edges: d-pbtn predicts regulations of igfr, mapk and cdk by erbb indeed, binding of erbb and igfr has been shown in lysate from skbr cells resistant to trastuzumab by immunoprecipitation. similarly, there is experimental evidence for an activation of mapk by erbb as we predicted. the inferred activation of cdk by erbb to our knowledge has not been reported before, and may warrant further experiments to confirm or falsify this prediction. some evidence also exists for the predicted strong activation of myc by esr. both proteins have been shown by affinity chromatography to bind to one another, and blocking of esr was recently shown to decrease myc expression. similarly, activation of akt increases ccnd and myc expression, and a regulation of ccne by myc and the activation of ccnd by mapk have been demonstrated experimentally. interestingly, d-pbtn predicts an activation of myc and ccne by igfr. both interactions to our knowledge have not been previously reported, and may be interesting candidates for experimental validation. discussion in this manuscript, we present d-pbtn, a novel method to infer signal-transduction networks from time course perturbation data. we evaluate d-pbtn in an extensive simulation study, and demonstrate its application on experimental data regarding erbb signaling. the simulation study shows stable inference results in light of missing data and noise, and reasonable robustness with respect to model hyperparameters. on the erbb application, we have run the inference with and without a network prior, and provide both results for comparison in figure while the overall topology is similar whether or not a prior network is used, there are significant differences in individual edges. the quality of network inference with a network prior clearly depends on the quality of the prior used as well as on the amount and quality of experimental data available, we therefore recommend to carefully evaluate the effect a given network prior has. in our example, there is experimental evidence for many of the predicted interactions, both from the run with and without network prior, directly confirming the quality of the inference. furthermore, in addition to confirming known interactions, d-pbtn predicts several novel edges, which may be interesting candidates for experimental validation. running time remains a major concern with d-pbtn, mainly due to the time required for sampling from the posterior distribution. each single sampling step requires running time o, where k is the number of perturbations, t is the number of time steps, and n is the number of nodes in the network. this can be easily seen from the likelihood function, which essentially iterates over all time points, perturbations and nodes in the network, and computes equation in each step, which again iterates over all nodes. running time for the likelihood dominates the time required for the prior distribution, hence the total running time for the posterior is o. unfortunately, the shape and complexity of the likelihood function, and thus the number of sampling steps required to sample from the posterior p, also depends on k, t and n in a nontrivial way, and the choice of required sampling steps is a nontrivial problem. clearly running time increases linearly with the number of sampling steps, but how exactly the minimum number of sampling steps required depends on k, t and n is unclear. while the present implementation of d-pbtn is a non-parallel implementation, evolutionary mcmc is straightforward to parallelize, allowing for a very efficient parallel implementation of the sampler. in practical terms, running time on the sn example was approximately hours with the present, non-parallelized version of the algorithm, and computations much beyond the size of the erbb network with its nodes appear computationally prohibitive without further parallelization. an efficient parallel implementation on a larger compute cluster in principle should make d-pbtn suitable also for "larger" network inference problems with several dozen, possibly up to a few hundred nodes, but this is clearly not an approach that is suitable for networks with thousands of proteins, let alone proteome-wide models. as with all bayesian approaches, prior distributions and prior hyperparameters impact results, and need to be carefully chosen. reasonable values for some of the parameters can be estimated from biological expectations, for example the expected level of "connectedness" in a network can be used to estimate parameters q and s of the prior on w, or information about the speed of signal transduction in a network could be used to estimate hyperparameters for the prior on τ. however, still a significant amount of "gut feeling" and experience are required to get good estimates. while our simulation study indicates reasonable stability with respect to parameter changes, a possible extension of our work is to employ empirical bayes’ approaches to set model hyperparameters, or, granted enough experimental data and compute time are available, crossvalidation-schemes could be used for hyperparameter estimation. the comparison between d-pbtn and the non time-series version pbtn on the different simulated data sets gave some very interesting results. while d-pbtn outperformed pbtn on all data sets when aucroc was used as performance measure, pbtn yielded higher median aucpr than d-pbtn on datasets with high levels of noise, or data simulated using a different underlying mathematical model. aucpr in contrast to aucroc is less influenced by an unequal distribution of class labels, the superior results of pbtn over d-pbtn under certain conditions are therefore of relevance for the choice of inference method. it can be shown that if a roc curve for method a dominates a roc curve for method b, then also the corresponding pr curve for method a dominates the pr curve for method b, and vice versa. the inverted ranking of d-pbtn and pbtn in aucroc and aucpr implies that the roc curves of d-pbtn and pbtn do not dominate one another. rather, there are certain domains in the roc and pr plots where one or the other method becomes better, and correspondingly the curves for d-pbtn and pbtn intersect. this in particular occured in our simulation study for the data with high levels of noise. it is tempting to speculate that under conditions with many "wrong" data points, focusing on a stable steady state may have advantages over trying to capture the full dynamics of the system, as d-pbtn does. this point will require further research in the future, to determine what exactly determines which of the two approaches is superior under given conditions, and to provide guidelines for method selection. on the other hand, if the interest is in feedback loops and cycles in the networks, pbtn should not be used, as this approach cannot infer cyclic networks. interestingly, on the erbb application example, d-pbtn clearly outperformed pbtn both with respect to aucroc and aucpr.  in this manuscript, we have developed a new method to infer signal transduction networks from time course perturbation data. based on relatively few time points after a large number of different perturbations, our method is able to reconstruct the underlying signaling network with high accuracy. the mathematical approach we employ is based on dynamical bayesian networks, and assumes discrete states and time steps, at which signaling molecules are either "on" or "off". the approach is therefore suitable to experiments where protein activation or protein expression is completely or almost completely switched off, such as with rnai knockdowns. several related approaches exist. our work extends pbtn by explicitly taking time into account. in contrast to pbtn, the method can thereby infer feedback and feedforward loops. provided sufficient time-resolved measurements after knockdowns are available, a unique network can be inferred using d-pbtn. this is in contrast to methods such as pbtn, depn and nested effects models, which only return equivalence classes or single networks out of an equivalence class. in this regard, our method is comparable to dynamic depns and dynamic nems, which also explicitly take time course data into account. however, d-nems have been developed for "effect" observations after knockdowns, which only provide very indirect information about the underlying network. d-depns, on the other hand, require long time series and work most efficiently if a relatively small number of perturbations are carried out. the typical experimental scenario encountered in practice is different, comprising short time courses after a large number of perturbations. d-pbtn have explicitly been developed for this situation. due to the bayesian approach pursued, prior biological knowledge can easily be integrated into d-pbtn inference, and marginalization over unobserved proteins or time points is straightforward and very efficient due to the likelihood simulation. we employ an mcmc approach to sample from the posterior distribution of models given the experimental data, which allows it to compute probability distributions over alternative network topologies. this is particularly useful if multiple models explain the data equally well, and can be used to design additional experiments to then distinguish further between high-scoring topologies. using halite clustering, similar models can be grouped together and probabilities estimated from the mcmc samples. overall, d-pbtn is a powerful approach for network inference when time-resolved measurements of the response of a dynamical system after perturbation are available. under these conditions, the approach outperforms other state-of-the-art methods. most approaches available in the literature consider either changes in steady state levels after perturbation, or the temporal dynamics of an unperturbed system for network inference. d-pbtn makes use of both the dynamical aspects and the perturbation effects in network reconstruction. we present an extensive simulation study to evaluate the approach, showing stable performance under varying conditions of noise and missing values in the data, and we demonstrate the application of d-pbtn to infer signaling networks in the erbb system. here, our approach both reconfirms known interactions, as well as suggesting some novel edges as candidates for further experimental study. availability and requirements project name: dynamic probabilistic threshold networks project home page:http://www.kaderali.org/software/dpbtn.html. operating system: tested on microsoft windows programming language: c# other requirements: none license: gnu gpl any restrictions to use by non-academics: none electronic supplementary material additional file 1: supplementary material. additional file one contains supplementary text with additional information. additional file 2: supplementary figure s effect of the amount of unobserved proteins on inference performance. shown is the distribution of the area under the precision recall curve for replicates of simulated data sets, over the fraction of unobserved proteins, for d-pbtn, depn, bdagl and pbtn on dataset sn the dashed line shows expected auc pr results for random guessing. additional file 3: supplementary figure s effect of noise on inference performance. the figure shows the distribution of auc pr values obtained for different levels of noise on sn for d-pbtn, depn, bdagl and pbtn. the dashed line shows expected auc pr results for random guessing. additional file 4: supplementary figure s inference performance on the reconstruction the kegg networks. this figure shows performance of the d-pbtn, depn bdagl and pbtnapproaches on subnetworks taken from the kegg database, with data simulated as described in methods. shown are auc pr values of simulated data sets, generated from ten different kegg subnetworks. the dashed line shows expected auc pr results for random guessing. additional file 5: supplementary figure s inference of negative feedback cycles: bean-plot showing the distributions of inferred edge weights. a unimodal distribution is observed with the full data, whereas no unique network can be identified when only observations of two downstream proteins are given. additional file 6: supplementary figure s effect of model parameters on performance: the figure shows the effect of varying the stochasticity parameter on the auc values. network inference has been done using data simulated from simulated network auc roc values are shown as red squares, auc pr values are shown as blue circles. additional file 7: supplementary figure s effect of model parameters on performance: the figure shows the effect of varying the rate and shape parameter q and s of the prior p on the auc values. network inference has been done using data simulated from simulated network panels a and b: effect of simultaneous changes of q and s by up to 20% on auc roc and auc pr. panels b and c: effect of changing parameters q and s individually by up to 50%. additional file 8: supplementary figure s effect of model parameters on performance: the figure shows the effect of varying the rate and shape parameters r and s of the prior on the bias weight w by up to 50%. network inference has been done using data simulated from simulated network auc roc values are shown as red squares, auc pr values are shown as blue circles. additional file 9: supplementary figure s experimental reference network extracted from string for erbb network inference. competing interests the authors declare that they have no competing interests. authors’ contributions lk and nak designed the research and methodology. nak and lk wrote the software and analyzed the data. lk and nak wrote the paper. both authors read and approved of the final version of the manuscript. we acknowledge funding from the german ministry of education and research via systec/fancy and gerontosys/agenet, and from the european union 7th framework program.
it exploits the dynamic response of a system after external perturbation for network reconstruction. on simulated data, we show that the approach outperforms current state of the art methods. we then use our method to study egf-mediated signaling in the erbb pathway. on the erbb data, our approach recovers a significant fraction of the known interactions, and predicts novel mechanisms in the erbb pathway. we utilize dynamic bayesian networks with probabilistic boolean threshold functions to describe protein activation. dynamic probabilistic threshold networks is a new method to infer signaling networks from time-series perturbation data. we present a novel method to infer signaling networks from time-course perturbation data. network inference deals with the reconstruction of molecular networks from experimental data. the model posterior distribution is analyzed using evolutionary mcmc sampling and subsequent clustering, resulting in probability distributions over alternative networks. due to data limitations, this typically is an ill-posed problem, and requires the integration of prior biological knowledge or strong regularization.
 network inference deals with the reconstruction of molecular networks from experimental data. given n molecular species, the challenge is to find the underlying network. due to data limitations, this typically is an ill-posed problem, and requires the integration of prior biological knowledge or strong regularization. we here focus on the situation when time-resolved measurements of a system’s response after systematic perturbations are available.  we present a novel method to infer signaling networks from time-course perturbation data. we utilize dynamic bayesian networks with probabilistic boolean threshold functions to describe protein activation. the model posterior distribution is analyzed using evolutionary mcmc sampling and subsequent clustering, resulting in probability distributions over alternative networks. we evaluate our method on simulated data, and study its performance with respect to data set size and levels of noise. we then use our method to study egf-mediated signaling in the erbb pathway.  dynamic probabilistic threshold networks is a new method to infer signaling networks from time-series perturbation data. it exploits the dynamic response of a system after external perturbation for network reconstruction. on simulated data, we show that the approach outperforms current state of the art methods. on the erbb data, our approach recovers a significant fraction of the known interactions, and predicts novel mechanisms in the erbb pathway. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. issue-copyright-statement© the author 2014 