Integrating prior knowledge in automatic network reconstruction

Abstract

The reconstruction of models from experimental data is a challenging problem due to the inherited complexity of biological systems. We developed an exact, exclusively data-driven approach to reconstruct Petri nets from experimental time-series data. Our approach aims at reconstructing all such netwoks that t the given experimental data, to provide all possible alternatives of mechanisms behind the experimental observations, which typically results in a large set of solution alternatives. To keep this solution set reasonably small while still guaranteeing its completeness, we rstly generate only Petri nets being minimal in the sense that all other networks tting the data contain the reconstructed ones. We further aim at avoiding the generation of minimal solutions which are \technically correct" but would be ruled out later during a subsequent veri cation process to check whether the returned solutions are \biological meaningful" or even contradict well-established biological knowledge. For that, we propose to extent the considered input (beyond the information given with the experimental time-series data) for the reconstruction process and demonstrate with the help of a running example the in uence on the generated solution set.

-

Systems biology aims at the integrated experimental and theoretical analysis of molecular or cellular networks to achieve a holistic understanding of biological systems and processes. To gain the required insight into the underlying biological processes, experiments are performed and experimental data are interpreted in terms of models. Depending on the biological aim, the type and quality of the available data, di erent types of mathematical models are used and corresponding reconstruction methods have been developed. Our work is dedicated to Petri nets which turned out to coherently model both static interactions in terms of networks and dynamic processes in terms of state changes, see e.g. [ 7,10 ].

In fact, a network P = (P; T; A; w) re ects the involved components by places p 2 P and their interactions by transitions t 2 T , linked by weighted directed arcs (p; t); (t; p) 2 A. Each place p 2 P can be marked with an integral number xp of tokens de ning a system state x 2 Zj+P j, i.e., we obtain X := fx 2 ZjP j : xp 0g as set of potential states. A transition t 2 T is enabled in a state x if xp w(p; t) for all p with (p; t) 2 A, we denote the set of all such transitions by T (x). Switching t 2 T (x) yields a successor state succ(x) = x0 with x0p = xp w(p; t) for all (p; t) 2 A and x0p = xp + w(t; p) for all (t; p) 2 A. Dynamic processes are represented by sequences of such state changes.

Our central question is to reconstruct models of this type from experimental time-series data by means of an exact, exclusively data-driven approach. This approach takes as input a set P of places and discrete time-series data X 0 X given by sequences (x0; x1; : : : ; xk) of experimentally observed system states. The goal is to determine all Petri nets (P; T; A; w) that are able to reproduce the data, i.e., that perform for each xj 2 X 0 the experimentally observed state change to xj+1 2 X 0 in a simulation. Hence, in contrast to the normally used stochastic simulation, we require that for states where at least two transitions are enabled, the decision between the alternatives is not taken randomly, but a speci c transition is selected. Thus, (standard) Petri nets have to be equipped with additional activation rules to force the switching of speci c transitions (to reach xj+1 from xj ), and to prevent all others from switching. For that, di erent types of additional activation rules are possible.

On the one hand, control-arcs are used to represent catalytic or inhibitory dependencies. An extended Petri net P = (P; T; (A [ AR [ AI ); w) is a Petri net which has, besides the (standard) arcs in A, two additional sets of so-called control-arcs: the set of read-arcs AR P T and the set of inhibitor-arcs AI P T ; we denote the set of all arcs by A = A [ AR [ AI . Here, a transition t 2 T (x) coupled with a read-arc (resp. an inhibitor-arc) to a place p 2 P can switch only if at least w(p; t) tokens (resp. less than w(p; t) tokens) are present in p; we denote by TA(x) the set of all such transitions.

On the other hand, in [ 9,12,13 ] priority relations among the transitions of a network are employed to re ect the rate of the corresponding reactions, where the fastest reaction has highest priority and, thus, is taken. In Marwan et al. [ 9 ] it is proposed to model such priorities with the help of partial orders O on the transitions. We call (P ; O) an (extended) Petri net with priorities, if P = (P; T; A; w) is an (extended) Petri net and O a priority relation on T . Priorities can prevent enabled transitions from switching: For each state x, a transition t 2 TA(x) is allowed to switch only if there is no other enabled transition in TA(x) with higher priority; we denote by TA;O(x) the set of all such transitions.

We call (P ; O) X 0-deterministic if TA;O(x) contains at most one element for each experimentally observed state x 2 X 0. Based on earlier results in [ 3,4,5,9,13 ], an integrative method to reconstruct all X 0-deterministic extended Petri nets with priorities tting given experimental time-series data X 0 was proposed in [ 6 ] (see Section 2).

Our approach aims at reconstructing all netwoks of the studied type that t the given experimental data, to provide all possible alternatives of mechanisms behind the experimentally observed phenomena. Typically, this results in a large set of solution alternatives. To keep this solution set reasonably small while still guaranteeing its completeness, we generate only Petri nets being minimal in the sense that all other networks tting the data contain the reconstructed ones. Here, we propose a method to insert only minimal sets of control-arcs during the reconstruction process (see Section 2). We further aim at avoiding the generation of minimal solutions which are \technically correct" but would be ruled out later during a subsequent veri cation process to check whether the returned solutions are \biological meaningful" or even contradict well-established biological knowledge. For that, we extent the considered input by integrating biological prior knowledge (beyond the information given with the experimental time-series data) into the reconstruction process and demonstrate with the help of a running example the in uence on the generated solution set (see Section 3). We close with some concluding remarks and lines of future work. 2

Reconstructing extended Petri nets with priorities We describe the input, the main ideas, and the output of our approach from [ 6 ]. Input. A set of components P (standing for proteins, enzymes, genes, receptors or their conformational states, later represented by the set of places) is chosen which is expected to be crucial for the studied phenomenon. If P contains known P -invariants (subsets P 0 P of places where the sum of the number of all tokens on all the places in P 0 is constant, e.g., di erent functional states of a cell or conformational states of a molecular complex), they are collected in a set IP .

To perform an experiment, one rst triggers the system in some state x0 (by external stimuli like the exposition to a pathogen), to generate an initial state x1. Then the system's response to the stimulation is observed and the resulting state changes are measured for all considered components at certain time points. This yields a sequence of (discrete or discretized) states xj 2 ZjP j re ecting the time-dependent response of the system to the stimulation in x1, which typically terminates in a terminal state xk where no further changes are observed. The corresponding experiment is X 0(x1; xk) = (x0; x1; : : : ; xk): Several experiments starting from di erent initial states in a set Xi0ni X 0, reporting the observed state changes, and ending at di erent terminal states in a set Xt0erm X 0 describe the studied phenomenon, and yield experimental time-series data of the form X 0 = fX 0(x1; xk) : x1 2 Xi0ni; xk 2 Xt0ermg: Thus, the input of the reconstruction approach is given by (P; IP ; X 0).

Example 1. As running example, we will consider experimental biological data from the light-induced sporulation of Physarum polycephalum as in [ 6,13 ]. In P. polycephalum plasmodia, the photoreceptor involved in the control of sporulation Spo is a protein which occurs in two stages PF R and PR. The developmental decision of starving P. polycephalum plasmodia to enter the sporulation pathway is controlled by environmental factors like visible light [ 11 ]. If the darkadapted form PF R absorbs far-red light F R, the receptor is converted into its red-absorbing form PR, which causes sporulation after several hours [ 8 ]. If PR is exposed to red light R, it is photo-converted back to the initial state PF R (photoreversion), which prevents sporulation if the red light pulse is given shortly after the far-red pulse, but not if the red pulse is delivered after more than an hour when the phytochrome photoreceptor has had su cient time to cause the formation of a biochemical downstream signal G that subsequently causes the sporulation of the cell. The changes between the stages PF R and PR only require fractions of seconds and can be experimentally observed due to a change of color of the phytochrome protein. The experimental setting consists of For a successful reconstruction, the data X 0 need to satisfy two properties: reproducibility and monotonicity. The data X 0 are reproducible if for each xj 2 X 0 there is a unique observed successor state succX 0 (xj ) = xj+1 2 X 0. Reproducibility is obviously necessary and can be ensured by preprocessing [ 15 ]. Whether a state xj 2 X 0 and its observed successor succX 0 (xj ) = xj+1 2 X 0 are also consecutive system states depends on the chosen time points to measure the states in X 0. In fact, xj+1 may be obtained from xj by a switching sequence of some length, where the intermediate states are not reported in X 0. The data X 0 are monotone if for each pair (xj ; xj+1) 2 X 0, the values of the elements do not oscillate in the possible intermediate states between xj and xj+1. It was shown in [ 4 ] that monotonicity has to be required, too (which is equivalent to demand that all essential responses are indeed reported in X 0). Output. An extended Petri net with priorities (P; O) with P = (P; T; A; w) ts the given data X 0 when it is able to perform every observed state change from xj 2 X 0 to xj+1 2 X 0. For that, associate with P an incidence matrix M 2 ZjP j jT j whose rows correspond to the places p 2 P and whose columns M t to the update vector rt of the transitions t 2 T : >8 w(p; t) if (p; t) 2 A; < rpt = Mpt := +w(t; p) if (t; p) 2 A;

>:0 otherwise: Reaching xj+1 from xj by a switching sequence using the transitions from a subset T 0 T is equivalent to obtain xj+1 from xj by adding the corresponding columns M t of M for all t 2 T 0, i.e., xj + Pt2T 0 M t = xj+1:

T has to contain enough transitions to perform all experimentally observed switching sequences. The underlying standard network P = (P; T; A; w) is conformal with X 0 if, for any two consecutive states xj ; xj+1 2 X 0, the linear equation system xj+1 xj = M has an integral solution 2 NjT j such that represents a sequence (t1; :::; tm) of transition switches, i.e., there are intermediate states xj = y1; y2; :::; ym+1 = xj+1 with yl +M tl = yl+1 for 1 l m. The extended Petri net P = (P; T; A; w) is catalytic conformal with X 0 if tl 2 TA(yl) for each intermediate state yl, and the extended Petri net with priorities (P; O) is X 0-deterministic if ftlg = TA;O(yl) holds for all yl.

The desired output consists of all minimal X 0-deterministic extended Petri nets (P; O) (all having the same set P of places as part of the input). Example 2. We represent in Fig. 3 (page 54) several alternative X 0-deterministic extended Petri nets tting the experimental data X 0 from our running example.

We now brie y sketch the reconstruction procedure.

Representing the observed responses. As initial step, extract the observed changes of states from the experimental data. For that, de ne the set

D := dj = xj+1

xj : xj+1 = succX 0 (xj ) 2 X 0 : Generating the complete list of all X 0-deterministic extended Petri nets P = (P; T; A; w) includes nding the corresponding standard networks and their incidence matrices M 2 ZjP j jT j. Due to monotonicity [ 4 ], it su ces to represent any dj 2 D using sign-compatible update vectors from the following set only: 8 > > Box(dj ) = < > > : r 2 ZjP j : 0 dp rp dp if djp > 0 >9 rp 0 if djp < 0 => rp = 0 if djp = 0 Pp2P 0 rp = 0 8P 0 2 IP ;>> n f0g: Next, we determine for any dj 2 D, the set (dj ) of all integral solutions of dj =

X rt2 Box(dj) trt; t 2 Z+; (1) and for each 2 (dj ), the (multi-)set R(dj ; ) = frt 2 Box(dj ) : t (6=dj0)garoef update vectors used for this solution . By construction, Box(dj ) and always non-empty since dj itself is always a solution due to reproducibility [ 6 ]. Every permutation = (rt1 ; : : : ; rtm ) of the elements of a set R(dj ; ) gives rise to a sequence of intermediate states xj = y1; y2; :::; ym; ym+1 = xj+1 with =

; (xj ; dj ) = (y1; rt1 ); (y2; rt2 ); : : : ; (ym; rtm ) :

To compose all possible standard networks, we have to select exactly one solution (dj ) for each dj 2 D and to take the union of the corresponding sets R(dj ; 2) in order to yield the columns M t = rt of an incidence matrix M of a conformal network. To ensure that the generated conformal networks can be made X 0-deterministic, we proceed as follows.

Detecting priority con icts between sequences. By construction, every sequence ; (xj ; dj ) induces a priority relation O , since it implies which transition ti is supposed to have highest priority (and thus switches) for every intermediate state yi. To impose valid priority relations O among all transitions of the reconstructed networks, we have to take con icts between priority relations O induced by di erent sequences into account. Two sequences and 0 are in priority con ict t = rt0 and intermediate states y; y0 such that if there are update vectors r 6 t; t0 2 T (y) \ T (y0) and (y; rt) 2 but (y0; rt0 ) 2 0 (since this implies t > t0 in O but t0 > t in O 0 ). We have a weak (resp. strong) priority con ict if y 6= y0 (resp. y = y0) which can (resp. cannot) be resolved by adding control-arcs.

We construct a priority con ict graph G = (VD [Vterm; ED [EW [ES ) whose nodes correspond to all possible sequences ; (xj ; dj ) and whose edges re ect weak and strong priority con icts (WPC and SPC for short): { VD contains the sequences ; (xj ; dj ) for all xj j = succX 0 (xj ) xj , for all 2 (dj ) and all permutat2ionXs0 n oXft0eRr m(dja;nd).d { Vterm contains for all xk 2 Xt0erm the trivial sequence (xk; 0). { ED contains all edges between two sequences ; 0 stemming from the same di erence vector. { ES (resp. EW ) re ects all SPCs (resp. WPCs) between sequences ; 0 stemming from distinct di erence vectors.

In G, we generate all node subsets S selecting exactly one sequence ; (xj ; dj ) per di erence vector dj 2 D such that no SPCs occur between the selected sequences and the nodes in Vterm. Clearly, a node 2 VD can never be selected for any solution S if it is in SPC with a terminal state sequence or with all sequences 0 stemming from one di erence vector dj 2 D. Removing all such nodes and their incident edges from G yields the reduced priority con ict graph G0. We can show that the sets of suitable selections S obtained in G and G0 are equal and that there is always at least one feasible selection.

Example 4. We obtain the following reduced priority con ict graph G0 for the running example, where bold edges indicate SPCs and thin edges WPCs.

Q3 σ1(x1,d1) σ3(x1,d1) σ1(x7,d3) σ1(x3,d3) σ(x8,0)

Q1 σ(x0,0) σ1(x6,d4) σ3(x6,d4) σ(x4,0) σ1(x5,d4) σ3(x5,d4) σ1(x2,d2)

Q3 From G0, we obtain as feasible subsets Si = S0 [ Si0 with

S0 = f 1(x2; d2); 1(x3; d3); 1(x7; d3)g S10 = f 1(x1; d1); 1(x5; d4); 1(x6; d4)g, S50 = f 3(x1; d1); 1(x5; d4); 1(x6; d4)g, S20 = f 1(x1; d1); 1(x5; d4); 3(x6; d4)g, S60 = f 3(x1; d1); 1(x5; d4); 3(x6; d4)g, S30 = f 1(x1; d1); 3(x5; d4); 1(x6; d4)g, S70 = f 3(x1; d1); 3(x5; d4); 1(x6; d4)g, S40 = f 1(x1; d1); 3(x5; d4); 3(x6; d4)g, S80 = f 3(x1; d1); 3(x5; d4); 3(x6; d4)g. Constructing X 0-deterministic Petri nets. Every selected subset S corresponds to a conformal standard network PS = (P; TS ; AS ; w): we obtain the columns of the incidence matrix MS of the network by taking the union of all sets R(dj ; ) corresponding to the sequences = ; (xj ; dj ) selected by 2 S. Example 5. We apply the method to the feasible sets S1 and S4 from Exp. 4 and obtain the standard networks presented in Fig. 2 with TS1 = D and in Fig. 3 with TS4 = fd1; d2; d3; r4:1; r4:2g, respectively.

If there are weak priority con icts 0 2 EW for nodes ; 0 2 S [ Vterm, denoted by WPC( ; 0), the constructed standard network PS needs to be made X 0-deterministic by inserting appropriate control-arcs. By [ 6 ], a WPC( ; 0) between two sequences and 0 involving update vectors rt 6= rt0 and intermediate states y 6= y0 with t; t0 2 T (y) \ T (y0) such that (y; rt) 2 but (y0; rt0 ) 2 0 has to be resolved by adding appropriate control-arcs that either turn rt into a transition t which is disabled at y0 (then t > t0 forces t to switch in y whereas only t0 is enabled at y0), or vice versa. Let P (y; y0) be the set of places where y and y0 di er and CA( ; 0) the set of all possible control-arcs that can resolve WPC( ; 0). Then CA( ; 0) partitions into two subsets CAt>t0 ( ; 0) disabling t at y0 containing { a read-arc (p; t) 2 AR with weight w(p; t) > yp0 8p 2 P (y; y0) with yp > y0 , p { an inhibitor-arc (p; t) 2 AI with w(p; t) < yp 8p 2 P (y; y0) with yp < yp0, and CAt<t0 ( ; 0) disabling t0 at y containing { a read-arc (p; t0) 2 AR with weight w(p; t0) > yp 8p 2 P (y; y0) with yp0 > yp, { an inhibitor-arc (p; t0) 2 AI with w(p; t0) < yp0 8p 2 P (y; y0) with yp0 < yp. Remark 1. If one of y; y0 is a terminal state, say y0, then CAt<t0 ( ; 0) = ; follows since t has to be disabled at y0 and t > t0 = 0 holds automatically. Moreover, if y = y0 then P (y; y0) = ; follows which is the reason why SPCs cannot be resolved by adding control-arcs.

Due to [ 6 ], inserting one control-arc from CA( ; 0) resolves the WPC( ; 0) in PS . Here, we further discuss mutual in uences of control-arcs in the resulting extended Petri nets as well as the issue of only constructing minimal catalytic conformal networks.

On the one hand, inserting a control-arc (p; t) 2 CA( ; 0) in PS might disable t at a state in another sequence 00 2 S n ; 0 where t is supposed to switch. In this case, (p; t) has to be removed from CA( ; 0), resulting in a reduced set CAS ( ; 0). On the other hand, one control-arc may resolve several WPCs in PS if the corresponding sets CAS ( ; 0) intersect.

Therefore, we propose the following consideration: Introduce one variable z(p;t) 2 f0; 1g for each possible control-arc (p; t) 2 CAS ( ; 0) for all WPCs in PS . Construct a 0=1-matrix AS whose columns correspond to all those variables (resp. control-arcs) and whose rows encode the incidence vectors of the sets CAS ( ; 0) for all WPCs in PS . Then any 0=1-solution z of the system AS z 1 encodes a suitable set of control-arcs resolving all WPCs in PS and, thus, a hitting set or cover of AS . By [ 14 ], we are only interested in nding minimal models tting X 0, where minimality is interpreted in the sense that all non-minimal models contain another one also tting the data. We can show that using nonminimal covers of AS yields non-minimal extended Peri nets but that we need all of them for the sake of completeness, which corresponds to determining the so-called blocker b(AS ) of AS .

Example 6. We list all WPCs between sequences in our running example: WPC1 between WPC2 between WPC3 between WPC4 between WPC5 between WPC6 between WPC7 between WPC8 between WPC9 between WPC10 between WPC11 between WPC12 between WPC13 between σ1(x2, d2) and σ1(x3, d3) and σ1(x7, d3) and σ3(x1, d1) and σ3(x1, d1) and σ3(x1, d1) and σ1(x2, d2) and σ3(x5, d4) and σ3(x5, d4) and σ3(x6, d4) and σ3(x6, d4) and σ1(x5, d4) and σ3(x5, d4) and σ(x0, 0) σ(x0, 0) σ(x0, 0) σ1(x7, d3) σ(x0, 0) σ(x8, 0) σ3(x5, d4) σ1(x3, d3) σ(x4, 0) σ1(x3, d3) σ(x4, 0) σ3(x6, d4) σ1(x6, d4) due to due to due to due to due to due to due to due to due to due to due to due to due to d2, 0 d3, 0 d3, 0 d3, r1.2 r1.2, 0 r1.2, 0 d2, r4.2 d3, r4.2 r4.2, 0 d3, r4.2 r4.2, 0 d4, r4.2 d4, r4.2 ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈

T (x0) ∩ T (x2) T (x0) ∩ T (x3) T (x0) ∩ T (x7) T (x1) ∩ T (x7) T (x0) ∩ T (x1) T (x8) ∩ T (x1) T (x2) ∩ T (x5) T (x3) ∩ T (x5) T (x4) ∩ T (x5) T (x3) ∩ T (x6) T (x4) ∩ T (x6) T (x5) ∩ T (x6) T (x5) ∩ T (x6) We obtain the following control-arcs to resolve WPCs between sequences: CA(WPC1) = CA(WPC2) = CA(WPC3) = CA(WPC4) = CA(WPC5) = CA(WPC6) = CA(WPC7) = CA(WPC8) = CA(WPC9) = CA(WPC10) = CA(WPC11) = CA(WPC12) = CA(WPC13) = {(PFR, d2) ∈ AI, (PR, d2) ∈ AR}, {(PFR, d3) ∈ AI, (PR, d3) ∈ AR, (G, d3) ∈ AR}, {(G, d3) ∈ AR}, {(F R, d3) ∈ AI, (G, d3) ∈ AR, (F R, r1.2) ∈ AR, (G, r1.2) ∈ AI}, {{((SFpRo,,rr11..22)) ∈∈ AARI,}(,G, r1.2) ∈ AI}, {(R, r4.2) ∈ AR, (R, d2) ∈ AI}, {(R, r4.2) ∈ AR, (G, r4.2) ∈ AI, (R, d3) ∈ AI, (G, d3) ∈ AR}, {(R, r4.2) ∈ AR, (Spo, r4.2) ∈ AI, (G, r4.2) ∈ AI}, {(R, d3) ∈ AI, (R, r4.2) ∈ AR}, {(R, r4.2) ∈ AR, (Spo, r4.2) ∈ AI}, {(G, d4) ∈ AI, (G, r4.2) ∈ AR}, {(G, d4) ∈ AR, (G, r4.2) ∈ AI}.

The following reductions of sets of possible control-arcs are necessary: for all standard networks PSi , we have 1(x3; d3) and 1(x7; d3) selected simultaneously, which both are in WPC with (x0; 0), see WPC2 and WPC3. To resolve WPC2, we have CA(WPC2) = f(PF R; d3) 2 AI ; (PR; d3) 2 AR; (G; d3) 2 ARg.

However, (PF R; d3) 2 AI and (PR; d3) 2 AR do not only disable d3 at x0, but also d3 at x7 (due to x0PF R = x7PF R = 1 and x0PR = x7PR = 0). Since d3 is supposed to switch at x7, we obtain CASi (WPC2) = f(G; d3) 2 ARg as reduced set of possible control-arcs to resolve WPC2 in all networks P(Si). Similarly, (G; r4:2) 2 AI has to be excluded from CA(WPC8) and CA(WPC9) in PS4 and PS8 as otherwise r4:2 would be disabled at x6 where it is supposed to switch by 3(x6; d4) 2 S4; S8.

For the feasible set S4, we obtain as matrix AS4 :

WPC1 WPC2 WPC3 WPC7 WPC8 WPC9 WPC10 WPC11 (PFR,dX2) ∈ AI (PR,d3X) ∈ AR (G,d3X)∈ AR (R,r4.2) ∈ AR (R,d2) ∈ AI (R,d3) ∈ AI (Spo,r4.2) ∈ AI

XX X

X X X X X

X X

X X The blocker b(AS4 ) contains four minimal covers of AS4 which correspond to the di erent sets of control-arcs in the four extended Petri nets depicted in Fig. 3, all arising from the standard network PS4 .

For S1, the matrix AS1 contains the rst 3 rows and columns of AS4 , the blocker b(AS1 ) contains two minimal covers of AS1 which correspond to the control-arcs in the two extended Petri nets in Fig. 2 arising from PS1 .

Note that b(AS ) is non-empty if and only if none of the sets CAS ( ; 0) is empty. We can show that there is at least one catalytic conformal network for any given X 0. All catalytic conformal extended Petri nets PS;P 0 = (P; TS ; AS;P 0 ; w) based on PS can be made X 0-deterministic by taking all the priorities O for all 2 S, where O is de ned by O = nti > t : t 2 TAS;P 0 (yi) n ti; 1 i mo : By construction, there are no priority con icts in the extended network PS;P 0 between O and O 0 for any ; 0 2 S, thus we obtain the studied partial order OS;P 0 = [ O : 2S R d4

PFR

d1

PR This nally implies: Theorem 1. Every extended network PS;P 0 = (P; TS ; AS;P 0 ; w) together with the partial order OS;P 0 is an X 0-deterministic extended Petri net, and there are no other minimal extended Petri nets with priorities tting the given data X 0. d4

PR

d1

PFR d2 d3 PS1b

G

Spo PFR d2 d3

PR PS4b

PR PFR d2 d3 PS4d r4.2 R r4.1 G Spo r4.2

G

Spo r4.1

FR d1

R PR

PFR d2 G d3 Spo Standard network PS4

FR d1 r4.2 R r4.1 r4.2 d1 FR

FR d1 r4.1

FR d1

Example 7. For the two extended Petri nets based on PS1 in Fig. 2, no priorities are needed to obtain X 0-deterministic extended Petri nets (PS1a ; ;) and (PS1b ; ;).

For the four extended Petri nets in Fig. 3 based on PS4 , the priority relation O4 = f(r4:2 > d2)g is required for PS4a and PS4b , whereas O4 = f(d2 > r4:2); (d3 > r4:2)g is required for PS4c and PS4d , to obtain X 0-deterministic extended Petri nets.

In total, the complete solution set contains 66 minimal X 0-deterministic extended Petri nets with priorities, see Table 1. 3.1

Integrating prior biological knowledge

Indecomposability of di erence vectors

In the step Representing the observed responses, the set of potential update vectors which might constitute the incidence matrices of the networks are considered. Hereby, for each dj 2 D, the set Box(dj ) contains only sign-compatible vectors due to monotonicity (which avoids homogeneous solutions in (1) according to minimality) and takes P -invariants into account (which avoids infeasible intermediate states according to prior biological knowledge). In some cases, one could restrict Box(dj ) further, e.g., if { dj exactly corresponds to a well-known biochemical reaction (including the correct stoichiometry) or to a well-known mechanism (that a certain trigger is detected by a suitable receptor), { experiments have shown that subsets of the input components of dj do not lead to the observed response, { dj is treated as black box-like reaction where only input and output matter, but not the intermediate mechanism due to the chosen level of detail. In such cases, we propose to exclude the corresponding response dj 2 D from decomposition and, instead, just de ne Box(dj ) := fdj g in accordance with the existing knowledge.

Example 8. Since the light-dependent reactions of the photoreceptor are so much faster than the subsequent processes that are considered in the reconstruction process, the di erence vector describing the phytochrome photoconversion will not be decomposed into di erent reaction vectors.

Thus, for the di erence vectors d1, d4 and d5 only the canonical sequences 1(x1; d1), 1(x5; d4) and 1(x6; d4) remain in the priority con ict graph, and S1 remains as only possible selection. Accordingly, the total number of solutions reduces from 66 to the 2 presented in Fig. 2, see Table 1. 3.2

Treating equal di erence vectors in the same way In the step Detecting priority con icts between sequences, all observed responses dj 2 D are treated independently from each other, so far. If, however, two di erence vectors di; dj 2 D are equal, we clearly have Box(di) = Box(dj ) and, thus, (di) = (dj ). Here, it is natural to require that both di; dj are decomposed in the same way (i.e., by the same 2 (di) = (dj )) and that the involved reactions are applied in the same order (i.e., by the same permutation of the elements in the sets R(di; ) = R(dj ; )) to obtain the same transitions (with equal control-arcs and priorities) in both cases. Indeed, using the same {

2 (di) = (dj ) corresponds to the fact that the same subset of molecules involved in a reaction will not interact according to di erent mechanisms, { permutation of the elements in R(di; ) = R(dj ; ) corresponds to the fact that the order in which the reactions are applied re ects the relative rates of the reactions in R(di; ) = R(dj ; ), so the same relation between reaction rates shall lead to the same priorities within the resulting sequences. We call two sequences ; (xi; di) and 0 ; (xj ; dj ) twin sequences if di = dj and the same and has been used. To force that twin sequences are always selected together, we propose to add strong priority con icts between all other sequences stemming from a pair di; dj of equal di erence vectors while creating the priority con ict graph, since no two sequences in strong priority con ict are selected for the same network.

Example 9. In our running example, we have d3 = d6 and d4 = d5. The latter vectors can be decomposed in di erent ways, among the resulting sequences we have 1(x5; d4), 1(x6; d5) and 3(x5; d4), 3(x6; d5) as pairs of twin sequences. Forcing to select these pairs together rules out the four selections S2; S3; S6; S7 so that only S1, S4, S5, S8 remain as possible selections. Accordingly, the total number of solutions reduces from 66 to 18, as reported in Table 1. 3.3

Knowledge on relative reaction rates

In the step Constructing X 0-deterministic Petri nets, to resolve a WPC( ; 0) between ; 0 involving update vectors rt 6= rt0 and intermediate states y 6= y0, either transition t has to be disabled at y0 or transition t0 at y, while the decision between t and t0 on the other state can be handled by a priority. Here, prior knowledge about the relative reaction rates of t and t0 (e.g. gained from the time-scales during the experiments) could help to decide whether t > t0 or t < t0 better re ects the reality, and than choose between control-arcs either from CAt>t0 ( ; 0) or from CAt0>t( ; 0).

So far, this idea is already applied for WPCs involving a terminal state: if y0 is a terminal state and 0 = (y0; 0) its trivial sequence, then t > 0 holds automatically at y, but t has to be disabled at y0 using control-arcs from CAt>0( ; 0) whereas CAt<0( ; 0) is empty.

This idea can be generalized to any WPC( ; 0) where the time-scale of the corresponding experimental observations clearly di ers in order to deduce the correct priority t > t0 or t < t0. In such cases, we propose to reduce CA( ; 0) accordingly either to CAt>t0 ( ; 0) or to CAt0>t( ; 0).

Example 10. In our running example, we have WPC1, WPC2, WPC3, WPC5, WPC6, WPC9, WPC11 involving a terminal state; the according reductions of the sets CA( ; 0) to CAt>0( ; 0) are already applied in Exp. 6.

Moreover, WPC4, WPC7, WPC8, WPC10 involve reactions with clearly different time-scales during the experimental observations: { d1 and d4 = d5 need only milliseconds to occur, { d2 needs about 1 hour to occur, and { d3 = d6 need at least 10 hours. Accordingly, we can reduce the sets CA( ; 0) as follows: { due to r1:2 > d3, for WPC4 only (F R; r1:2) 2 AR and (G; r1:2) 2 AI remain; { due to r4:2 > d2, for WPC7 only (R; r4:2) 2 AR remains; { due to r4:2 > d3, for WPC8 only (R; r4:2) 2 AR and (G; r4:2) 2 AI remain, while for WPC10 only (R; r4:2) 2 AR is left.

At least one of these WPCs occurs in the standard networks coming from the selected sets S2 S4; S6 S8. Note that in the two remaining WPCs the timescale of the involved responses is equal (see Exp. 6). Consequently, the number of extended Petri nets decreases from 66 to 36 as reported in Table 1. 4

Discussion The subject of this paper was an approach from [ 6 ] that aims at reconstructing all X 0-deterministic extended Petri nets that t given experimental data X 0, to provide all possible alternatives of mechanisms behind the experimentally observed phenomena. This typically results in a large set of solution alternatives. To keep this solution set reasonably small while still guaranteeing its completeness, we rstly generate only Petri nets being minimal in the sense that all other networks tting the data contain the reconstructed ones. In the presented approach, the minimality concept is applied twice: { monotone data: using only sign-compatible vectors in Box(dj ) avoids homogeneous solutions during the decomposition of (dj ) (and super ous transitions in the networks), see [ 4 ]. { minimal hitting sets : we here proposed to use only minimal sets of controlarcs to resolve all weak priority con icts in a standard network which avoids unnecessary control-arcs (and arti cial dependencies), see Section 2. This ensures that the presented approach exactly generates all minimal extended Petri nets with priorities (Theorem 1). We further avoid generating minimal solutions which are \technically correct" but would be ruled out later during a subsequent veri cation process to check whether the returned solutions are \biological meaningful" or even contradict well-established biological knowledge as in [ 2 ]. For that, we extent the considered input by integrating prior biological knowledge beyond the information given with the experimental time-series data into the reconstruction process. We propose to integrate prior knowledge in the following way: { P-invariants : helps to obtain feasible intermediate states in all sequences which avoids the generation of solutions contradicting known facts, see [ 13 ]. { indecomposable di erence vectors : help to keep already known subnetworks or mechanisms, see Section 3.1. { treating equal di erence vectors in the same way : helps to keep consistency in the interpretation of experimental observations, see Section 3.2. { obeying terminal states and reaction rates : helps to chose meaningful priorities and to avoid arti cial control-arcs, see Section 3.3.

So far, the algorithmic procedure due to [ 6 ] takes as input (P; IP ; X 0) where the list of P -invariants and the monotonicity of X 0 are already used in the rst step Decomposing di erence vectors to determine Box(dj ). To integrate further knowledge in the reconstruction procedure to force the algorithm to make \the right decisions" in some intermediate steps, we suggest, based on the previous discussions, to extend the input to (P; IP ; X 0; Din; Deq; OD) where { Din contains all indecomposable di erence vectors (for that, carefully select them according to the before mentioned critera, e.g., to preserve known subnetworks or mechanisms or to take a certain level of detail into account); { Deq contains all pairs of equal di erence vectors that shall be treated in the same way (for that, only chose equal di erence vectors where also the time elapsed during the experimental observation was equal); { OD contains information about the reaction rates between di erence vectors (for that, only impose priorities for su ciently di erent rates according to clearly di erent time scales during the experiments).

Example 11. The three examples from Section 3 can be interpreted as follows: { Exp. 8 shows the result taking (P; IP ; X 0; Din = fd1; d4; d5g; ;; ;) as input; { Exp. 9 shows the result with (P; IP ; X 0; ;; Deq = fd3 = d6; d4 = d5g; ;); { Exp. 10 the result with (P; IP ; X 0; ;; ;; OD = fd1; d4; d5 > d2 > d3; d6g). Whereas the rst setting reduces the number of solution alternatives to 2, the combination of the two latter scenarios reduces the number of solution alternatives to 12, see Table 1 below.

minimal Din 6= ∅ Deq 6= ∅

OD 6= ∅ Deq, OD 6= ∅

To conclude, we notice that providing indecomposable di erence vectors has the largest impact on the solution set. However, even if no indecomposable difference vectors can be identi ed, treating equal di erence vectors in the same way and deducing relative reaction rates from the time-scale of the experimental observations leads to a substantial reduction of the solution set, keeping only \biological meaningful" network alternatives.

The further goal is to provide an implementation for the presented reconstruction method, including the option of integrating prior knowledge as additional input during the reconstruction. For that, we will use Answer Set Programming as done for the reconstruction of standard networks in [ 1 ]. Finally, we plan to apply the presented reconstruction approach to di erent biological experimental data. We expect an important impact of Automatic Network Reconstruction on the integrated experimental and theoretical analysis of biological systems towards their holistic understanding, since computational models derived from experimental data by our exact, exclusively data-driven approach have predictive ability due to completeness of the solution set guaranteed by mathematical proofs.

References

1. M. Durzinsky , W. Marwan , M. Ostrowski , T. Schaub , and A. K. Wagler . Automatic network reconstruction using ASP . Th. and Pract . of Logic Progr., 11 : 749 { 66 , 2011 .

2. M. Durzinsky , W. Marwan , and A. K. Wagler . Reconstruction of extended Petri nets from time series data and its application to signal transduction and to gene regulatory networks . BMC Systems Biology, 5 , 2011 .

3. M. Durzinsky , W. Marwan , and A. K. Wagler . Reconstruction of extended Petri nets from time-series data by using logical control functions . Journal of Mathematical Biology , 66 : 203 { 223 , 2013 .

4. M. Durzinsky , A. K. Wagler , and R. Weismantel . A combinatorial approach to reconstruct Petri nets from experimental data . Lecture Notes in Comp. Sc. (special Issue CMSP 2008 ), 5307 : 328 { 346 , 2008 .

5. M. Durzinsky , A. K. Wagler , and R. Weismantel . An algorithmic framework for network reconstruction . J.Theor. Computer Science , 412 ( 26 ): 2800 { 2815 , 2011 .

6. M. Favre and A. K. Wagler . Reconstructing X 0-deterministic extended Petri nets from experimental time-series data X 0 (extended abstract) . CEUR Workshop Proceedings (Special Issue BioPPN 2013 ), 988 : 45 { 59 , 2013 .

7. I. Koch and M. Heiner . Petri nets . In B. H. Junker , F. Schreiber, editor, Biological Network Analysis , Wiley Book Series on Bioinformatics , pages 139 { 179 , 2007 .

8. T. Lamparter and W. Marwan . Spectroscopic detection of a phytochrome-like photoreceptor in the myxomycete physarum polycephalum and the kinetic mechanism for the photocontrol of sporulation . Photochem Photobiol , 73 : 697 { 702 , 2001 .

9. W. Marwan , A. K. Wagler , and R. Weismantel . A mathematical approach to solve the network reconstruction problem . Math. Methods of Operations Research , 67 ( 1 ): 117 { 132 , 2008 .

10. J. W. Pinney , R. D. Westhead , and G. A. McConkey . Petri net representations in systems biology . Biochem. Soc. Tarns. , 31 : 1513 { 1515 , 2003 .

11. C. Starostzik and W. Marwan . Functional mapping of the branched signal transduction pathway that controls sporulation in Physarum polycephalum . Photochem. Photobiol., 62 : 930 { 933 , 1995 .

12. L. M. Torres and A. K. Wagler . Encoding the dynamics of deterministic systems . Math. Methods of Operations Research , 73 : 281 { 300 , 2011 .

13. A. K. Wagler . Prediction of network structure . In I. Koch, F. Schreiber , and W. Reisig, editors, Modeling in Systems Biology , volume 16 of Computational Biology, pages 309 { 338 , 2010 .

14. A. K. Wagler and J.-T. Wegener . On minimality and equivalence of Petri nets (extended abstract) . Fundamenta Informaticae , 128 : 209 { 222 , 2013 .

15. A. K. Wagler and J.-T. Wegener . Preprocessing for network reconstruction: Feasibility test and handling infeasibilty (extended abstract) . CEUR Workshop Proceedings (Special Issue CS &P 2013 ), ( 1032 ): 434 { 47 , 2013 .