Representing dynamic biological networks with multi-scale probabilistic models

Dynamic models analyzing gene regulation and metabolism face challenges when adapted to modeling signal transduction networks. During signal transduction, molecular reactions and mechanisms occur in different spatial and temporal frames and involve feedbacks. This impedes the straight-forward use of methods based on Boolean networks, Bayesian approaches, and differential equations. We propose a new approach, ProbRules, that combines probabilities and logical rules to represent the dynamics of a system across multiple scales. We demonstrate that ProbRules models can represent various network motifs of biological systems. As an example of a comprehensive model of signal transduction, we provide a Wnt network that shows remarkable robustness under a range of phenotypical and pathological conditions. Its simulation allows the clarification of controversially discussed molecular mechanisms of Wnt signaling by predicting wet-lab measurements. ProbRules provides an avenue in current computational modeling by enabling systems biologists to integrate vast amounts of available data on different scales.

In addition we validated the knockdown efficiency on protein level (right panels) by Western-Blot (A-B) as well as by mass spectrometry (C-D). GAPDH (glycerinaldheyd-3-phosphat-dehydrogenase) was used as loading control for PCRs and Western-Blots. In the case of JNK1, the protein amount of JNK1 was below the detection limit.  Figure  1A is marked by a black box. GAPDH (glycerinaldehy-3-phosphat-dehydrogenase) was used as a house keeping gene. Note that samples were run on different geles due to the size of APC. The stealthRNA (stRNA) used is depicted above the images.

Introduction to ProbRules
A ProbRules model consists of a graph of interactions and a set of rules. The vertices of the graph represent components of a system. In models of biochemical systems, these can correspond to macromolecules like proteins, DNA and RNA or small molecule compounds like guanosine diphosphase (GDP). Possible interactions among these components are represented by edges in the graph. For a thorough introduction to graph theory see [1]. Besides direct interactions, such edges can also represent component states. For example, we used edges to represent protein modifications like phosphorylation in our study. Fig. N1 contains an example of a static interaction graph. To represent a state S t at time t of the model, probabilities p t are assigned to the interactions [2]. These probabilities represent a measure of the presence of the corresponding interactions. Accordingly, the presence of an interaction can range from 0, representing complete absence, to 1, representing full presence. These states S t can be seen as random graph models in which each edge is sampled with the associated probability. Furthermore, as the probabilities of the interactions are assigned to distinct time points t, their sequences can represent dynamic behavior.
To represent the dynamics of the biochemical network in the model, interdependencies between the interactions are encoded by rules. Each rule affects a single target interaction driving its state towards a defined target probability by a specified attack rate and the activity of the rule. The rule's activity at time t is evaluated as the probability that its Boolean formula holds on state S t using the rules of probability calculus. The interaction's state is also driven to its initial probability via a global decay rate when the rules for that interaction are not active. The example from Fig. N1 can be extended by interdependency rules as shown in Fig. N2. It includes two rules. In each of these rules a single source interaction triggers activation respectively inhibition of a target interaction. Now follows a more formal description of the ProbRules modeling approach.

Static Interaction Graph
The graph of interactions GI = (V, E) consists of a set V of vertices representing system components and a set E ⊆ V × V of edges representing possible interactions of the components. As such, the graph represents the static structure of the modeled system. The example shown in Fig

Model states, time points and dynamics
The state of an interaction (i, j) ∈ E of two components i ∈ V and j ∈ V at a time point t is denoted by the probability p t (i, j). A state S t (E) of the model for a time point t is defined by corresponding probabilities p t attached to the edges E of the graph GI: Each such S t defines a random graph model which essentially is a probability distribution D t over possible subgraphs G = (V, E G ) of GI with E G ⊆ E [3]. Therefore the probability P r(G) of a subgraph The edges (i, j) can be viewed as independent random variables that are true with probability p t (i, j). Dynamics of a modeled system can be represented by a sequence of states S 0 , S 1 , ..., S T . Thereby, the probabilities of the different interactions can evolve over time. Although the random variables corresponding to the edges are independent from each other at any particular point in time t, interdependencies between the interactions at different time points can be introduced. This is achieved in a controlled way by evaluating rules.

Interdependency Rules
A set R of rules defines the interdependencies between activity states of interactions. Each rule takes the form where φ represents a Boolean condition (formula) on source interaction states and (i, j) is the affected target interaction whose activity state p(i, j) is driven towards the target probability q by the attack rate a. Interdependencies between interactions can act by driving target interactions' probabilities towards arbitrary values using arbitrary attack rates. The generality of this formulation enables a straightforward translation of biochemical experimental findings into a ProbRules model as this was experienced during the development of the Wnt signaling pathway model. For the present work it was sufficient to use a single conjunction (Boolean AND) of the states of several, possibly negated source interactions in the formula φ.
For the extended example in Fig. N2 two rules can be specified: Given a model state S t at time t, the probability of activation φ t of a rule's formula φ can be determined using the rules of probability calculus. Thus, φ t denotes the probability that the logical formula φ holds in a randomly sampled subgraph according to the distribution D t .
The interdependency rules operate in a step-wise manner in regard to the evolution in time. At each point in time t − 1, each rule r proposes a new probability value q t (r, i, j) for its target interaction (i, j). First, a discussion of the effects of rules independently of one another follows and then a description how the effects of different rules are combined.

Interaction states proposed by rules
A basic rule r with condition φ, target (i, j), target probability q, and attack rate a has to be read as: whenever φ holds at time t − 1 the probability of the interaction (i, j) at time t will be That means that whenever φ t−1 = 1 the new probability for (i, j) is an affine combination of the previous state p t−1 (i, j) and the target probability q as determined by the attack rate a.
In general (for arbitrary S t ) the condition φ will not hold on all subgraphs that can be sampled according to D t , and thus the probability φ t will be < 1. To account for this, the rules contribute only with a corresponding factor φ t−1 · q t (r, i, j) to the target p t (i, j). Due to subgraphs where φ t is false, a factor for the negation (1 − φ t−1 ) would be lost at each time step. To account for this lost factor and subgraphs, a standard decay rule is introduced.

Decay
The standard decay rule states that an interaction that is not affected by any active rule returns to its defined initial state p 0 (i, j) by a global decay rate d:

Combining one rule with the decay rule
Combining a single effective rule with the default decay rule yields the value Thus the decay rule applies only when the other rule is not active.

Combining an arbitrary number of rules
Consider now rules r 1 , ..., r n with formulas φ 1 , ..., φ n which target the same interaction (i, j). They each propose new states q t (r 1 , i, j), ..., q t (r n , i, j) of that interaction. A method for combination can be obtained by considering the meaning of a rule r k . Such a rule basically states that whenever φ k holds -for those subgraphs for which φ k holds -the target has to be set to q t (r k , i, j). A concrete realization consists of subgraphs for which exactly m conditions φ k1 , ..., φ km hold with corresponding q t (r k1 , i, j), ..., q t (r km , i, j). In these cases, the combination will be the average to these subgraphs: As for the single rule case, the decay rule still applies for subgraphs in which n i=1 φ i is false.

Synchronous update
After deriving the proposed next states p t of all interactions (i, j) based on the previous interaction state S t−1 yields the new interaction state S t . Then, a new cycle can be started which yields new interaction states and so on, until some final time point T is reached. This allows to simulate dynamics of biochemical systems like the Wnt signaling networks using their static interaction graph and interdependency rules on the states of the interactions.

Inputs and perturbations
Inputs can be provided to a ProbRules model at specific interactions and times by specifying an explicit probability This also allows to investigate perturbations like inhibition and constitutive activation of a specific interaction (i, j) in a specified ProbRules model.

Grounding on Causal Probabilistic Time Logic
The introduced ProbRules model is based on notions from probabilistic logic programming [4] and statistical relational learning [5]. In particular, it is inspired on the CPT-L framework which forms a logical foundation. This section describes how the ProbRules interdependency rules (like CPT-L) can be translated into ProbLog [3], a recently introduced probabilistic extension of the programming language Prolog [6].

Causal Probabilistic Time Logic (CPT-L)
The ProbRules interdependency rules are a variation on the rules used in the CPT-L formalism of [7]. 1 CPT-L rules are of the form where the h i are propositional (Boolean) variables, φ is a Boolean expression over these variables, and p i ∈ [0, 1] probabilities with n i=1 p i = 1. The semantics of a CPT-L rule is that whenever the expression φ is true at time t − 1, one of the h i will be sampled with probability p i and will become true at time t. As we are working in a probabilistic logic the expression φ will only be true with a particular probability p φ . This probability will then be taken into account in the sampling process: the h i 's will be sampled with probability p φ · p i .
A CPT-L theory is a set of CPT-L rules. It may have multiple rules targeting the same propositional variables. These rules will be combined as two independent alternative causes to generate these targets. This corresponds to a disjunctive interpretation of the rules.

Stochastic Relational Process
A CPT-L theory specifies a probability distribution over sequences of logical interpretations. An interpretation is an assignment of truth-values to all propositional variables in the theory. So, a CPT-L theory defines a probability distribution over sequences of such interpretations I 1 , ..., I t which is a stochastic relational process.

Semantics of interdependency rules
The ProbRules interdependency rules defined before are essentially variants of CPT-L rules. Summarizing on this: • The propositions represent the presence or absence of an edge in the interaction graph. This implies that an interpretation corresponds to a subgraph of the interaction graph.
• The target probabilities of a CPT-L theory are immediate, whereas those of the ProbRules interdependency rules may involve more than one time step. Nevertheless, in the previous section was shown how the immediate probabilities can be computed.
• Each ProbRules interdependency rule has only a single target, whereas a CPT-L theory may have multiple targets (each h i is a target).
• Different ProbRules interdependency rules are combined using averaging whereas CPT-L rules are combined using OR.

Implementing ProbRules interdependency rules
Among the ProbRules interdependency rules that can easily be mapped onto a set of CPT-L rules are the cases with only a single rule for a target (i, j): which can be represented in CPT-L by Here each CPT-L rule has only one possible conclusion (p(i, j)). The probabilistic constants a, d and q are represented as Boolean variables that are true with their corresponding probability. Evaluating the probability p(i, j) using these CPT-L rules corresponds to computing the probability of the logical formula: Since the four terms are disjoint, the probability of the disjunctive expression is the sum of the probabilities of the expressions. The resulting expression has the same value as in Equation 1.
A set of n ProbRules interdependency rules for the same target (i, j) can be represented in a similar way in CPT-L. The two rules for decay are retained yielding as before: where Φ = ∨ i φ i . The attack rules are now of the following form: ..,km} are as defined in Equation 3. There is a pair of such rules for each possible ki in any possible subset {k1, ..., km} ⊆ {1, ..., n}. The m ki are defined using the following rule: which assures that the m ki are mutually exclusive and each of them has the correct probability 1/m. For any fixed {k1, ..., km} exist 2m rules for p(i, j). Each pair of rules will contribute 1 m q t (r ki , i, j) to p t (i, j) when Φ {k1,...,km} holds. This yields the same result as in Equation 2. Because the formulas Φ {k1,...,km} are mutually exclusive, the q t (Φ {k1,...,km} , i, j) can simply be summed to account for the attacking rules. This is combined with the decay rule.
This basically shows that a set of ProbRules interdependency rules inherits its function from CPT-L. It defines a probability distribution over finite sequences of subgraphs. As already indicated, particular interest is in the sequence of states (S t (E)) t≥0 that is produced using a ProbRules model (GI = (V, E), R), its rule set R, its interaction graph GI = (V, E) and given an initial state S 0 . Using CPT-L, each S t = {p t (i, j)|(i, j) ∈ E} is represented as a set of probabilistic propositions. The next state is then determined by computing the probability of (i, j) at time t + 1 using the above CPT-L rules. In this way, the sequence of states S t imposes a stochastic graph process over discrete time points.

Implementation in Prolog
Effective techniques for probabilistic inference have been developed and integrated in the programming language Prolog [6,3]. The ProbRules implementation requires SWI-Prolog, which is available for all major operating systems. Alternatively, one can also use the SWI-Prolog Docker container available from Docker Hub. Implementations of several network motifs and the Wnt model are available at Github.
Further details are available in the README file there. See also SWI-Prolog Quickstart for details on using SWI-Prolog.

Specifying ProbRules models
ProbRules models are specified based on Prolog. For a more elaborate introduction to Prolog see [8,6].

How to use the ProbRules implementation for model simulation
In order to use the ProbRules implementation for the simulation of a model of interest, we need to define involved interactions, interdependency rules between these and parameters which are used. At the start, we load the ProbRules implementation from ProbRules.pl: Parameters for target probabilities and attack rates are defined directly. The ProbRules implementation uses the global_decay parameter, and therefore it should be defined using that name. Other parameter can be defined as required for the rules:      Interactions consist of ordered pairs of components. The attached probability is used for initialization and as the target value for decay:   Inputs by interactions can be defined using the fixed predicate. We will use this in the Wnt model. Thus, fixed can be used for perturbations like knockout and overexpression of compounds, inhibition and constitutive activation of interactions in ProbRules models of biological systems.
In this example, we define a single simple regulation rule for the target interaction using no source, defined attack and target probability facts and a label which is required for distinguishing the rules:

Empty brakets [] enable a constantly active rule.
For autoregulation motifs we provide both activation and autoregulation rules:  For executing the simulation up to the time point 100, we use evaluation: 18 :-evaluation(100). 19 :-halt.
halt ends the Prolog interpreter and thus allows analyses to be scripted. Omitting halt opens an interactive environment that can be used for investigation of the model. The numbered lines can be stored in a file, e.g. srnarpar.pl. See also the interaction graphs in the next chapter. Running without any preambles or frills is achieved by typing swipl srnarpar.pl at a shell console or terminal. Timings were obtained on a recent MacBook Pro (Model 2018, i9, 2.9 GHz).   :-discontiguous (::)/2.

Network motifs interaction graphs and implementations
The following simulations reproduce the behavior of different network motifs introduced in [9,10]. The investigations included simple regulation, negative and positive auto-regulation, two variants of the bi-fan motif, coherent and incoherent feed-forward loops and the single input module. We represent activated forms of a component A by A* in the interaction graphs and by "as" in input files.
In the ODE implementations, regulation was carried out using Hill functions: Hill(X, K, n) = X n X n + K n It runs from 0 to 1 when input is provided as X and from 1 to 0 for input at K. X = K results in Hill(X, K, n) = 1 2 . The parameter n determines the slope of the sigmoidal curve.             :-evaluation(100).

Single input module
W* Figure N9: The single input module is a simple pattern in which one regulator activates a group of target genes with different strengths and thus speeds.

Using immediate rates disables ProbRules models to reproduce network motifs dynamics
In order to investigate advantages of the ProbRules models of network motifs over logical circuits we simplified the original network motifs implementations circumventing intermediate interactation activation probabilities. Instead, new target states are assigned immediately. This corresponds to setting all attack and decay rate values to 1. Such simplifications for our ProbRules models of network motifs resulted in dynamics depicted in Fig. N10. Here, we specified the same interactions and interdependency rules as in the ProbRules implementations used for Fig. 2 and changed all attack and decay rates to 1.
The symmetric bifan and the incoherent type I feed-forward loop fail to show any activation at the corresponding outputs. Simple regulation, negative and positive autoregulation as well as the coherent type I feed-forward loop and the single input module fail to represent differential dynamics between the outputs. Remarkably, in the asymmetric bifan both effects can be observed for the outputs resulting in dynamics contradicting both ODE and ProbRules implementations with intermediate rates.
In summary, only by using intermediate interaction activations the ProbRules models are able to represent network motifs.  Fig. 2. ProbRules models of network motifs on right side with the same interactions and rules as in Fig. 2, but with immediate rates.

List of abbreviations
The following logical relations and rules describing Wnt/β-catenin and Wnt/JNK signaling use these abbreviations:

In absence of a Wnt ligand
Logical relation 1 In absence of Wnt, Fz does not interact with LRP.
Speed fast Comment Both Fz [11,12,13,14,15] and LRP-5/6 [16,17,18,19] have been shown to interact with Wnt and to be responsible for Wnt induced signal transduction. Fz and LRP do not interact in absence of a Wnt ligand [20]. It has been shown that Wnt interacts with both, Fz and LRP that are brought in close proximity by this binding [15,18,20,21].

Logical relation 2 In absence of Wnt, Fz interacts with GαGDP and Gβγ
Speed average Comment Fz receptors are seven-pass transmembrane receptors. A whole series of publications indicate that Wnt signaling is likely G-protein coupled [22,23,24,25,26,27,28]. Data from [26] strongly suggest a direct interaction of Fz receptors with heterotrimeric G-proteins. Logical relation 2 represents the common view of signaling through heterotrimeric G-protein coupled receptors.

Logical relation 3
In absence of Wnt, LRP does not interact with axin.

Speed average
Comment LRP has been shown to interact with axin upon Wnt stimulation [29,30]. Wnt proteins induce the translocation of axin to the plasma membrane and increase the interaction of axin with LRP. This logical relation formalizes the situation in absence of a Wnt ligand.

Logical relation 4
In absence of Wnt, LRP does not interact with CKIγ or GSK3.

Speed fast
Comment LRP has been shown to interact with CKIγ and GSK-3β [31,32]. Both enzymes are able to modify LRP. This logical relation formalizes the situation in absence of a Wnt ligand.

Logical relation 5
In absence of Wnt, dsh does not interact with Fz.

Speed fast
Comment Dsh can interact with Fz [33,34,35]. It is thought, that this interaction occurs after Wnt stimulation, i.e. dsh is translocated to the membrane upon ligand binding [36,37] (see below for further details).

Logical relation 6
In absence of Wnt, the destruction complex forms in the cytosol and interacts with β-catenin.

Speed fast Comment
These steps describe the formation of the mature destruction complex, considering that GSK-3β also phosphorylates APC, which is favorable for the degradation of β-catenin [44]. This set of rules also indicates that the destruction complex is stable in case of APC modification. Axin also can be phosphorylated by GSK-3β [47,48] and this phosphorylation was shown to be required for the interaction between axin and β-catenin [49].

Speed fast b) If CKIα interacts with β-catenin, then β-catenin is modified to form
β-cat*. Speed fast c) In the presence of the mature destruction complex, GSK3 can interact with β-cat* and modify it into β-cat** (that is GSK phosphorylates β-cat* a total of three times).

Speed average
Comment This logical relation describes that β-catenin is marked for degradation through ubiquitination and subsequently degraded in the proteasome [38,51,52]. As a pre-requisite for ubiquitination, β-catenin is phosphorylated through CKIα and GSK-3β in a dual kinase mechanism [36,50]. Initially, β-catenin is phosphorylated by CKIα at Ser45. This modification is represented in this model as β-cat*. After this pre-phosphorylation, GSK-3β modifies residues Ser33, Ser37, and Thr41. This modified β-catenin is represented by β-cat**. This phosphorylated variant of β-catenin can be ubiquitinated through action of E1/2/3 to form β-catUb, which then is degraded through the proteasome.

Logical relation 9
If β-catenin is accumulating and interacts with JNK2, β-catenin can be translocated in the nucleus and interact with LEF.

Speed slow Comment
This logical relation indicates that free β-catenin which is not degraded through the proteasome can interact with transcription factors of the TCF/LEF family [53,54,55]. Moreover, recently it has been shown that RAC and JNK activity are also required for canonical wnt signaling activity [56,57,58], which is in accordance with our own experimental findings. This cooperative effect likely involves RAC1 (and maybe RAC3) and JNK2 and may be mediated via control of β-catenin stabilization [59] or nuclear localization [58]. Another possibility is a potential function of c-jun [57] or RAC (together with its Rho-GEF Tiam1) [56] as transcriptional co-activators of TCF/LEF.
Remark: This positive effect of RAC and JNK on canonical signaling activity is also expressed by logical relation j19 b).

Logical relation 10
If β-catenin in the nucleus interacts with LEF, LEF interacts with the DNA.

Speed average
Comment This logical relation formalizes the observation, that β-catenin functions as a transcriptional co-activator for TCF/LEF proteins [55,60,61,62,63,64,65]. The interaction of LEF/β-catenin with DNA results in a transcriptional response that was used as a read out for the pathway.
In presence of Wnt  This logical relation indicates the formation of an activated form of dishevelled, named dsh*. It has been shown that this requires activated Gα subunits [26]. The precise mechanism of this activation is unknown so far but likely involves CKIγ. CKIγ is one kinase that phosphorylates dsh and is required for Wnt induced stabilization of β-catenin [66,67]. Dishevelled has also been shown to interact with Gβγ [68]. If inactivated GαGDP interacts again with Gβγ before activation of dsh, no downstream signaling can occur.

Hint
Logical relations 12/13 and 14-16 are active at the same time.

Comment
For transduction of the Wnt signal through the Wnt/β-catenin signaling pathway, Wnt interacts with both, LRP and Fz. This logical relation indicates that the different branches through LRP and Fz occur simultaneously.

Speed average
Inhibition of phosphorylation b) If dsh* is formed, dsh* can interact with GBP/FRAT. Speed average c) If dsh* interacts with axin and GBP, GBP interacts with GSK3.
Speed average d) If GSK3 interacts with GBP, then GSK3 does not interact with axin, APC or β-catenin.

Speed fast
Comment There are two branches described, how β-catenin phosphorylation in the destruction complex can be inhibited. This set of rules summarizes how β-catenin phosphorylation is inhibited through inactivation of GSK-3β via dsh and GBP/FRAT. Dsh is required for Wnt signal transduction [36,40,66,69]. The phosphorylated form of dsh (likely phosphorylated by CKIγ) has a high affinity to GBP/FRAT that can interact with GSK-3β and can inactivate this enzyme [67,70,71]. The consequence of GSK-3β inactivation is that β-catenin accumulates; logical relations 9 and 10 will become effective!

---or ---Inhibition of ubiquitination
Remark: This logical relation does not make use of logical subrelations 17 b-d.

Comment
The signaling mechanism through inhibition of β-catenin ubiquitination as postulated by Li et al. [50] does not assume inhibition of β-catenin phosphorylation through inactivation of GSK-3β.
Logical relation 18 a) If dsh* is formed, dsh* can interact with Fz.
Speed average c) If dsh* interacts with Fz and axin, and LRP*interacts with Fz and axin, then axin interacts with PP1. Speed average Comment Phosphorylated axin can interact with LRP, in particular with the activated form LRP* [30,31,32,37,49]. Dsh interacts with Fz receptors [21,33] and is required for proper translocation of axin to the membrane [37].

Speed fast
Inhibition of phosphorylation a) If the destruction complex is inhibited and β-catenin is not degraded, β-catenin can accumulate.

Speed average
This set of rules summarizes the second possibility how β-catenin phosphorylation is inhibited, through a weakening of the interaction between β-catenin and axin, when axin is dephosphorylated by PP1 in response to the Wnt signal [49].

Speed average
Remark: Upon association of axin with phosphorylated LRP at the membrane, dissociation of E1/2/3 from the destruction complex was observed in one recent study [50]. Upon this dissociation of E1/2/3 from the destruction complex, β-catenin is no longer ubiquitinated and degraded [72].

Consequence of both versions of logical relations 18d) and 19a) is that
β-catenin will accumulate and therefore logical relations 9 and 10 will become effective. Concerning logical relation 19a) knockdown data obtained in this study hints, that Rac1 plays a role in β-catenin stabilization (independent of JNK) as a loss of Rac1 (but not JNK) diminishes the increase in free β-catenin after activation of canonical Wnt signaling (see also logical relation j17). This influence could be mediated e.g. via the guanine exchange factor DOCK4 [59].
Remark: Logical relation 19a) in both versions now assumes a positive influence of Rac on the amount of free β-catenin that is necessary for the model to fit the experimental data obtained in this study.
b) If β-catenin accumulates, the high β-catenin concentration will feed back to re-activate the destruction complex.

Speed average
Comment Two recent studies [73] as well as our own experimental data show, that upon Wnt stimulation the amount of phosphorylated β-catenin is markedly decreased for a short time but subsequently is recovered quickly. Based on simulation results obtained with an ODE model of β-catenin synthesis, phosphorylation and degradation, Hernandez et al. [74] hypothesize, that this transient drop occurs due to incomplete inactivation of the destruction complex, which thus can be re-activated at high β-catenin levels.
If one of the scaffold proteins axin or APC is missing, no destruction complex can be formed.

Speed fast Comment
This logical relation formalizes that the scaffold proteins axin and APC are critical components of the destruction complex. They are necessary to bring GSK-3β and β-catenin in close proximity thus enabling β-catenin phosphorylation (reviewed e.g. in [73]).
Logical relation 21 Independent of the presence or absence of Wnt, β-catenin is newly synthesized.

Speed average Comment
To be able to model the comparably fast turnover of β-catenin appropriately, besides the degradation through the proteasome (see logical relation 8f), also the β-catenin synthesis is modeled explicitly as in previous ODE-based models of Wnt signaling (e.g. [74,75]).

In absence of a Wnt ligand
Logical relation j1 In absence of Wnt, dsh does not interact with Fz.

Speed fast
Comment Dsh can interact with Fz [21,33,35]. It is thought that this interaction occurs upon Wnt stimulation, i.e. dsh is translocated to the membrane upon ligand binding [33,35].

Logical relation j2
In absence of Wnt, Fz does not interact with Ror2.

Speed fast Comment
Ror2 was found to associate with Fz and participate in JNK activation upon Wnt stimulation [76].

Logical relation j3
In absence of Wnt, Fz interacts with GαGDP and Gβγ.

Speed average
Comment A series of publications [22,23,25,26,27] indicate that Wnt signaling is likely G-protein coupled. Data from [26] strongly suggest a direct interaction of Fz receptors with heterotrimeric G-proteins.

Logical relation j4
In absence of Wnt, part of the Rac-GTPases is already present in the active, GTP-bound form. Speed average Comment It is known that small RhoGTPases of the Rho and Rac family and the subsequent MAPK (mitogen-activated protein kinase) cascade can be activated by different cytokines, coupling to tyrosine-kinase and G-protein coupled receptors (see e.g. [77,78]). Therefore it has to be assumed that even in the absence of Wnt ligands a certain amount of Rac will be present in the active GTP-bound state as long as other cytokines are present (e.g. from FCS (fetal calf serum) in the cell culture medium).

In presence of Wnt
Logical relation j5 In the presence of Wnt, Wnt interacts with Fz and Ror2.

Speed average
Comment It has been shown that Wnt interacts with both Fz [15] and Ror2 [76], with the latter two forming a complex [76].
Logical relation j6 a) If Wnt interacts with Fz, GαGDP is transferred into GαGTP. Gβγ still interacts with Gα at this step. Speed fast Comment GαGTP can be inactivated to GαGDP. GαGDP can interact again with Gβγ and Fz. This logical relation summarizes a common mechanism how signaling through heterotrimerc G-proteins can be turned off again.

Logical relation j8
If active GαGTP is released, GαGTP interacts with dsh.

Speed average
Comment It has been shown that the activation of dsh requires activated Gα subunits [26], although the precise mechanism for this activation is unknown so far. If inactivated GαGDP interacts again with Gβγ before activation of dsh, no downstream signaling can occur.

Logical relation j9
If G dissociates into GαGTP and Gβγ and Gβγ does not interact anymore with Gα or Fz, Gβγ interacts with PLCβ.

Speed average
Comment Gβγ has been shown to activate PLCβ [79]. If Gβγ interacts again with GαGDP before activation of PLCβ no downstream signaling can occur [79].

Logical relation j10
If free Gβγ interacts with and activates PLCβ, PLCβ forms DAG.
Logical relation j11 a) DAG can interact with PKCδ.
Speed average b) If DAG interacts with PKCδ, PKCδ is activated to form PKCδ*.

Speed fast
Comment As described in [81], PKCδ belongs to the PKC subfamily of novel PKCs. This class of PKCs generally binds to and is activated by DAG (comp. [81]) that is produced on the membrane upon extracellular signaling.

Logical relation j12
If PKCδ is activated by DAG to form PKCδ*, PKCδ can interact with dsh.

Speed average
Comment PKCδ was found to form a complex with dsh [81] and dsh also was shown to be phosphorylated by PKCδ [82]. Logical relation j12 reflects the essential role PKCδ plays in the Wnt/JNK pathway by activating dsh [81].
Logical relation j13 If Wnt interacts with Fz and Ror2, Ror2 interacts with Fz.
Speed average Comment Ror2 has been shown to interact with Wnt-5a [76,83,84] and to form a complex with Fz [76]. Also Wnt is known to interact with Fz (see logical relation j5) [15].

Logical relation j14
If Fz interacts with Ror2, Fz interacts with dsh.
Logical relation j15 If dsh interacts with Fz and dsh interacts with GαGTP and dsh interacts with PKCδ, dsh is activated to form dsh**. Speed average Comment It has been shown that the activation of dsh requires activated Gα subunits (see also logical relation j8) [26]. Moreover PKCδ -which forms a complex with dsh -also has been found to be important for the signaling from Fz to dsh [81].

Logical relation j16
If dsh** is formed, dsh can interact with Rac.
Speed average Comment Dsh has been observed to form a Wnt-induced complex with Rac [87,88] which has been found to function downstream of dsh in planar polarity signaling [89].
Logical relation j17 a) If dsh interacts with Rac, Rac is activated to form RacGTP.
Speed fast b) Only active GTP-bound Rac can interact with JNK1 and JNK2 and influence β-catenin.

Speed average
Remark: logical relations j17 b) and c) include a connection from Rac to β-catenin (possibly via DOCK4, but not including JNK) that is necessary for the model to fit the experimental data obtained in this study.

Comment
Logical relations j17 a) and b) summarize the common view how small RhoGTPases (like Rac) are activated: They are transferred from an inactive GDP-bound to an active GTP-bound state by Rho-GEFs in response to a stimulus and this GTP-bound form is required for further signal transduction. Concerning logical relation j17b) and c) knockdown data obtained in this study hints, that Rac1 plays a role in β-catenin stabilization independent of JNK, as a loss of Rac1 (but not JNK) diminishes the increase in free β-catenin after activation of canonical Wnt signaling. Therefore a connection between Rac and β-catenin independent of JNK has to be assumed in the model for the simulation results to fit the experimental data. This link could e.g. function via the guanine exchange factor DOCK4 [59].

Logical relation j18
If Rac interacts with JNK1 and JNK2, they are activated to form JNK1* and JNK2*.

Speed fast Comment
Logical relations j17 and j18 indicate that dsh, more exactly the c-terminal half of dsh which contains the DEP (EGL-10 and pleckstrin) domain [85], has been found to activate the JNK pathway [85,86] with the small GTPase Rac acting downstream of dsh. The signal transduction from dsh to JNK, reflected in logical relations j17 and j18, likely occurs through Rac, as a full correlation between JNK and Rac activation has been observed in different studies [77,78,87]. Possibly MLK3 serves as a link between Rac and JNK as shown in [90]. However there also exists data considering Rac unlikely to play a significant role in dsh-induced JNK activation [69].
Logical relation j19 a) If JNK1 and JNK2 are activated, they can interact with β-catenin.
Speed average b) If JNK2 interacts with free β-catenin, β-catenin can interact with LEF.

Comment
It has been shown that Rac and JNK activity are required for canonical Wnt signaling activity [56,57,58] which is in accordance with our own experimental findings. This cooperative effect likely involves Rac1 (and maybe Rac3) and JNK2 and may be mediated via control of β-catenin stabilization [59] or nuclear localization [58]. Another possibility is a potential function of c-jun [57] or Rac (together with its Rho-GEF Tiam1) [56] as transcriptional co-activators of TCF/LEF.
Remark: This positive effect of RAC and JNK on canonical signaling activity is also expressed in logical relation 10.
Speed average d) If JNK1 interacts with β-catenin and GSK3, and β-catenin is in the destruction complex, β-catenin can be phosphorylated by GSK3. Speed average Comment Contrary to the above mentioned findings there also exists data that activation of JNK can inhibit the transcription of canonical Wnt target genes by preventing β-catenin accumulation through induction either of β-catenin nuclear export [91] or β-catenin phosphorylation by GSK-3β [92].

Wnt model feedbacks
The ProbRules model of Wnt signaling contains a range of feedbacks. Fig. N11A contains all 69 interactions of the model represented as boxed rectangles. Each interaction in the Boolean formula of a rule targeting an interaction was depicted by an arrow. We found several subgraphs consisting of interactions that are all connected via loops to each other. The largest 5 are highlighted by colors. Fig. N11B

Reaching a stable state before Wnt stimulation
Our ProbRules model of Wnt signaling includes an initial maturation phase of the destruction complex. In this phase, the complex is not yet fully stable and thus also not yet fully active. During this initial phase of complex stabilization the output of the Wnt model shows an initial peak until the destruction complex is fully active and the output reaches a stable state after about 75 time steps. This is demonstrated in Fig. N12. However, this initial stabilization phase cannot be captured with experimental measurements in the used cell culture system where we always have to assume the presence of some fully mature destruction complexes in the absence of Wnt stimulation. In order to start our simulations of Wnt signal transduction from these defined conditions, we supplied the Wnt stimulus in our model only after 100 time steps. Additionally, for all further analyses we will concentrate on the network behavior after this initial stabilization phase and provide the simulation results starting from time step 76.  Figure N12: The output shows a small initial peak before Wnt stimulation due to a not fully established destruction complex.

Effect of the decay speed of the input signal
The input signal to the Wnt signaling network model was specified by predetermining the probabilities of Wnt ligand edges. The stimulus effectively consisted of a square pulse followed by an exponential decay: with t 0 = 100, t k = 175 and s = 0.02. In order to investigate the influence of the decay speed (s) of the input signal, we ran the simulation using a range of values for the corresponding parameter, i.e. the input decay rate (in the model implementation, the corresponding parameter is called slow_attack) . Fig. N13 shows the dynamics of the input interaction depending on the decay speed. Fig. N14 shows the resulting behavior of the output interaction (LEF/β-catenin -DNA). Although low values of slow_attack lead to a decrease of the decay speed at the input, for moderate and higher values of slow_attack we did not see a qualitative difference.

Effect of the β-catenin synthesis rate
In a further analysis, we investigated the influence of the β-catenin synthesis rate on the output signal. Therefore, we ran the simulation using a range of values for the corresponding parameter (called syn_rate in the model implementation). Fig. N15 shows the resulting behavior of the output interaction (LEF/β-catenin -DNA). Here, turning off β-catenin synthesis (syn_rate= 0) leads to a lower response. The response is not completely abolished as a certain basal amount of β-catenin is assumed to be present initially in our model to mimic the experimental situation in our cell culture system. There, a certain basal amount of β-catenin is also present initially [74,75]. As the behavior at the output does not change over a wide range (0.03..

Effects of additional rules
In this investigation, the original ProbRules model of Wnt signaling was perturbed systematically by adding artificial rules. Each of these rules is activated by one of the existing interactions in the network or a constitutively present interaction intended for permanent activation. They act on each of the non-input interactions in the network by either driving them towards full presence or full absence. This procedure results in 9380 new models, which were simulated and compared to the original specification.   Figure N16: Effects of a complete exploration by adding rules. Besides the original behavior, three major types of aberrant behavior were observed, each of which is shown in the plot by means of two representative curves. Additionally, for each type the percentage of additional rules inducing the respective behavior is given.
As mentioned in the main text, 89% of the 9830 added rules did not induce significant changes in the transcriptional Wnt response, i.e. after addition of each of these 8358 rules, the output of the modified models showed basically the same behavior as in the original model. Each of the remaining 11% of additional rules, on the other hand, induced some change in the transcriptional Wnt response. According to the elicited change in the transcriptional response these 1022 rules can be roughly divided in three types (see also main text Fig. 6D, reproduced here as Fig. N16), which we want to discuss in the following: Type 1: Rules of this type induce a reduction of the transcriptional Wnt response. Of the 576 rules of this type, 503 only reduce the strength of the signal (Type 1A), while 53 almost completely abolish the output signal (Type 1B). The observed reduction or elimination of the output signal can be caused through the inhibition of a target edge that is required for signal transduction or the activation of a target edge that inhibits signal transduction. For example, a reduction of the signal was caused by the addition of a rule through which the presence of the Wnt ligand induces the interaction of GSK3-β with β-catenin (example "Type 1A" in Fig. N16). An almost complete elimination of the output signal was for example caused by the addition of a rule through which the presence of the Wnt ligand inhibits the phosphorylation of dsh (example "Type 1B" in Fig. N16).
Type 2: Rules of this type induce an increase of the signal strength. An increase of the signal strength can be caused for example by introducing a shortcut through which an upstream source edge directly activates the output edge or by introducing a positive feedback loop through which an upstream target edge is activated by a downstream source edge. The introduced shortcut or feedback loop then leads to a stronger up-regulation and/or slower downregulation of the output signal. Of the 259 rules of this type, 171 cause only a moderate increase of the signal strength (Type 2A), while 88 rules induce a strongly increased output signal (Type 2B). A moderate increase (Type 2A) was caused for example by addition of a rule through which the interaction between TCF/LEF and the DNA feeds back to induce the interaction between JNK2 and β-catenin (example "Type 2A" in Fig. N16). A strongly increased output signal (Type 2B), on the other hand, was caused for example by addition of a rule through which the presence of the Wnt ligand directly induces the interaction between β-catenin and TCF/LEF (example "Type "2A" in Fig. N16).

Type 3:
Rules of this type induce TCF-dependent transcriptional activity independent of the Wnt stimulation. This Wnt-independent transcriptional activity can occur for example, if a source interaction that is present in the absence of Wnt activates a target interaction that in turn can activate downstream signaling. In contrast to rules of Type 2, in this case, the Wnt signal is not required for initiation of the response. Depending on the source and target edges, the strength of the observed effect can vary considerably. Such a Wnt-independent transcriptional activity was caused for example by addition of a rule through which the interaction between axin and CK1α induces the phosphorylation of dsh (example "Type 3A" in Fig. N16) or a rule through which an interaction between axin and β-catenin induces the presence of free β-catenin (example "Type 3B" in Fig. N16).
Strikingly, of the 67 interactions that occur as target edges in the additional rules (i.e. all network interactions except the two input signals), only 20 were observed to be involved in rules that caused abnormal network behavior. Rules in which the remaining 47 edges appeared as target edges did not affect the network behavior, irrespective of the source edge involved. Such sensitive target interactions were, besides the LEF-DNA interaction and the free β-catenin, mainly components of the destruction complex, the Wnt-receptor complex and interactions involving the heterotrimeric G proteins.