This work aims at reconstructing Petri net models for biological systems from experimental time-series data X 0. The reconstructed models shall reproduce the experimentally observed dynamic behavior in a simulation. For that, we consider Petri nets with priority relations among the transitions and control-arcs, to obtain additional activation rules for transitions to control the dynamic behavior. The contribution of this paper is to present an integrative reconstruction method, taking both concepts, priority relations and control-arcs, into account. Our approach is based on previous works for special cases and shows how these known steps have to be modi ed and combined to generate the desired integrative models, called X 0-deterministic extended Petri nets.
The overall aim of systems biology is to analyze biological systems and to
understand di erent phenomena therein as, e.g., responses of cells to
environmental changes, host-pathogen interactions, or e ects of gene defects. To gain the
required insight into the underlying biological processes, experiments are
performed and the resulting experimental data are interpreted in terms of models.
Depending on the biological aim and the type and quality of the available data,
di erent types of mathematical models are used and corresponding methods for
their reconstruction have been developed. Our work is dedicated to Petri nets,
a framework which turned out to coherently model static interactions in terms
of networks and dynamic processes in terms of state changes, see e.g. [
Our central question is to reconstruct models of this type from
experimental time-series data by means of an exact, exclusively data-driven approach.
We base our method on earlier results from [
On the one hand, in [
On the other hand, in [
For consistently integrating both concepts, priority relations and control-arcs,
into the modeling framework, the di culty is that both are concurrent concepts
to force or prevent the switching of enabled transitions. In [
A standard or simple Petri net P = (P; T; A; w) is a weighted directed bipartite graph with two kinds of nodes, places and transitions. The places p 2 P represent the system components (e.g. proteins, enzymes, genes, receptors or their conformational states) and the transitions t 2 T stand for their interactions (e.g., chemical reactions, activations or causal dependencies). The arcs in A (P T ) [ (T P ) link places and transitions, and the arc weights w : A ! N re ect stoichiometric coe cients of the corresponding reactions.
Each place p 2 P can be marked with an integral number xp of tokens, and any marking de nes a state x 2 NjP j of the system. In biological systems, all components can be considered to be bounded, as the value xp of any state refers to the concentration of the studied component p 2 P , which can only increase up to a certain maximum cap(p). This leads to a capacitated Petri net (P; cap), i.e., a Petri net P = (P; T; A; w) together with a capacity function cap : P ! N, whose set of potential states is X := fx 2 NjP j j xp cap(p)g: A transition t 2 T is enabled in a state x 2 X of a capacitated Petri net if E1 xp w(p; t) for all p with (p; t) 2 A, and, E2 xp + w(t; p) cap(p) for all p with (t; p) 2 A and we de ne T (x) := ft 2 T : t satis es E1, E2 in xg: 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 control-arcs by AC = AR [ AI , and the set of all arcs by A = A [ AR [ AI .
In a capacitated extended Petri net, switching of transitions is additionally controlled by read- and inhibitor-arcs; a transition t satisfying E1 and E2 can switch only if also the following conditions hold: E3 xp w(p; t) for all p with (p; t) 2 AR, and, E4 xp < w(p; t) for all p with (p; t) 2 AI .
In an extended Petri net, a transition is enabled in a state x 2 X if it satis es E1, : : : , E4 (otherwise, it is disabled ). The switch of a transition t enabled in x leads to a successor state succX (x) = x0 2 X whose marking is obtained by >8xp < >:xp;
w(p; t);
x0p :=
xp + w(t; p);
for all p with (p; t) 2 A;
for all p with (t; p) 2 A;
otherwise:
In general, there can be more than one transition satisfying E1, : : : , E4 in
a state x 2 X and we de ne TA(x) := ft 2 T : t satis es E1, : : : , E4 in xg:
The decision which transition switches is typically taken randomly (and the
dynamic behavior is analyzed in terms of reachability, starting from a certain initial
state). This is not appropriate for modeling biological systems which show a
deterministic behavior, e.g., where a certain stimulation always results in the same
response. In this case, additional activation rules are required in order to force
the switch from a state x to a speci c successor state succX (x). For this purpose,
priorities between the transitions of the network can be used to determine which
of the transitions in TA(x) has to be taken. Note that these priorities typically
re ect the rate of the corresponding reactions where the fastest reaction has
highest priority. In Marwan et al. [
Note that priorities can prevent enabled transitions from switching: for a state x 2 X , only a transition t 2 TA(x) is allowed to switch or can switch if E5 there is no other transition t0 2 TA(x) with (t t0) 2 O.
The set of all transitions that are allowed to switch in x is denoted by
TA;O(x) := ft 2 T : t satis es E1, : : : , E5 in xg:
To enforce a deterministic behavior, TA;O(x) must contain at most one element
for each x 2 X to enforce that x has a unique successor succX (x), see [
In this paper we consider capacitated extended Petri nets with priorities (P; cap; O): extended Petri nets P = (P; T; A; w) with a capacity function cap : P ! N on their places and a partial order O T T on their transitions. Our goal is to reconstruct X 0-deterministic extended Petri nets from given experimental data X 0. 3
Reconstructing X 0-deterministic extended Petri nets In this section, we describe the input, the main ideas, and the generated output of our integrative reconstruction approach. 3.1
A set of components P (later represented by the set of places) is chosen which is expected to be crucial for the studied phenomenon. All known P -invariants 1 of the system (e.g., di erent conformational stages of a cell, a receptor, a protein) shall be collected in a set IP .
To perform an experiment, one rst triggeres the system in some state x0 (by external stimuli like the change of nutrient concentrations or the exposition to some pathogens), to generate an initial state x1. Then the system's response to the stimulation is observed and the resulting state changes are measured 1 Laxly said, a P-invariant is a set P 0 P of places (components) where the sum of the number of all tokens on all the places in P 0 is constant. P-invariants are not computed by the algorithm but must be known a priori by a biologist. for all 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 for all components p 2 P at certain time points, 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. The
developmental decision of starving P. polycephalum plasmodia to exit the vegetative
plasmodial stage and to enter the sporulation pathway is controlled by environmental
factors like visible light [
P = fF R; R; PF R; PR; Spog IP = fPF R; PRg
X 0(x1; x3) = (x0; x1; x2; x3) X 0(x4; x0) = (x2; x4; x0)
XtXerinmi == ffxx13;; xx40gg
0 0 as input for the algorithm, we represent all observed states schematically in Fig 1.
x xFR xR xxPPFRR xSpo
x0 0 FR 0 1 00 d4
x1 1 0 1 0
0 0 1 0 1 0 x4 d1 R
x2 0 0 0 1 0 d2
x3 0 0 0 1 1
In the best case, two consecutively measured states xj ; xj+1 2 X 0 are also consecutive system states, i.e., xj+1 can be obtained from xj by switching a single transition in T . This is, however, in general not the case (and depends on the chosen time points to measure the states in X 0), but xj+1 is obtained from xj by a switching sequence of some length, where the intermediate states are not reported in X 0.
For a successful reconstruction approach, 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. Moreover, the data X 0 are monotone if for each pair (xj ; xj+1) 2 X 0, the possible intermediate states xj = y1; y2; :::; ym+1 = xj+1 satisfy yp1 yp1 yp2 yp2 : : : : : : ym p ym p ym+1 for all p 2 P with xj p p ypm+1 for all p 2 P with xjp xj+1 and
p xjp+1.
Whereas reproducibility is obviously necessary, it was shown in [
A capacitated extended Petri net with priorities (P; cap; 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 succX 0 (xj ) = xj+1 2 X 0. This can be interpreted as follows. With P, an incidence matrix M (P) 2 ZjP j jT j is associated, where each row corresponds to a place p 2 P of the network, and each column M (P) t to the update vector rt of a transition t 2 T : rpt = M (P)pt := >8 w(p; t) <
+w(t; p) >:0 if (p; t) 2 A; if (t; p) 2 A; otherwise: Reaching xj+1 from xj by a switching sequence using the transitions from a subset T 0 T is equivalent to obtain the state vector xj+1 from xj by adding the corresponding columns M (P) t of M (P) for all t 2 T 0: xj + X M (P) t = xj+1:
t2T 0 Hence, T has to contain enough transitions to perform all experimentally observed switching sequences. The underlying standard network P = (P; T; A; w) is 2 This is equivalent to say that system states in X 0 have been measured to appropriate time points such that the values of their components do not oscillate between two measured states or, equivalently, that all essential responses are indeed reported in the experiments. conformal with X 0 if, for any two consecutive states xj ; succX 0 (xj ) = xj+1 2 X 0, the linear equation system has an integral solution 2 NjT j such that is the incidence vector of a sequence (t1; :::; tm) of transition switches, i.e., there are intermediate states xj = y1; y2; :::; ym+1 = xj+1 with yl + M (P) 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 of the reconstruction approach consists of the set of all X 0-deterministic extended Petri nets (P; cap; O) (all having the same set P of places and the same capacities cap deduced from X 0). 3.3
To solve the problem of representing the observed responses by switching
sequences, we propose the following approach, based on previous works in [
D := dj = xj+1
xj : xj+1 = succX 0 (xj ) 2 X 0 : Example 2. From our running example in Fig. 1 we obtain D = fd1; d2; d4g with d1 = x2 x1 = ( 1; 0; 1; 1; 0)T , d2 = x3 x2 = (0; 0; 0; 0; 1)T and d4 = x0 x4 = (0; 1; 1; 1; 0)T .
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.
The rst step is to describe the set of potential update vectors which might constitute the columns of M .
Representation of di erence vectors. Recall that two consecutively measured states xj ; xj+1 2 X 0 are not necessarily consecutive system states, i.e., xj+1 may be obtained from xj by a switching sequence of some length, where the intermediate states are not reported in X 0. Due to monotonicity, the values of the elements cannot oscillate in the intermediate states between xj and xj+1.
Moreover, for any P -invariant P 0 2 IP , all suitable update vectors have to
satisfy Pp2P 0 rp = 0. Hence, it su ces to represent any dj 2 D using only
vectors from the following set
8
>
>
Box(dj ) = <
>
>
:
r 2 ZjP j :
0
djp
rp djp if djp > 0 9>
rp 0 if djp < 0 =>
rp = 0 if djp = 0
Pp2P 0 rp = 0 8P 0 2 IP ;>>
n f0g:
Remark 1. In previous approaches [
Next, we determine for any dj 2 D, the set (dj ) of all integral solutions of the equation system
X rt2 Box(dj) dj =
trt; t 2 Z+:
By construction, Box(dj ) and (dj ) are always non-empty since dj itself is always a solution due to the required reproducibility of the input data X 0 (which particularly includes dj 6= 0 for all dj 2 D). For each 2 (dj ), construct the (multi-)set
R(dj ; ) = frt 2 Box(dj ) : t 6= 0g of update vectors used for this solution .
Example 3. For D from Exp. 2, the update vectors for a decomposition are Box(d1) = fd1; r1; r2g, Box(d2) = fd2g and Box(d4) = fd4; r3; r4g with vectors r1 = ( 1; 0; 0; 0; 0)T , r2 = (0; 0; 1; 1; 0)T , r3 = (0; 1; 0; 0; 0)T and r4 = (0; 0; 1; 1; 0)T . Hence, the possible decomposition of the responses are d1 = d1 = r1 + r2, d2 = d2 and d4 = d4 = r3 + r4 and the resulting sets are R(d1; 1) = fd1g; R(d1; 2) = fr1; r2g; R(d2; ) = fd g
2 ;
R(d4; 1) = fd4g; R(d4; 2) = fr3; r4g: 3.4
To compose all possible standard networks, we have to select exactly one solution 2 (dj ) for each dj 2 D and to take the union of the corresponding sets R(dj ; ) 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.
Sequences and their con icts. 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 ) : By construction, every such sequence respects monotonicity and 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 priority con icts between priority relations O induced by di erent sequences into account.
Two sequences and 0 are in priority con ict if there are update vectors rt 6= rt0 and intermediate states y; y0 such that 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 priority con ict (WPC) if y 6= y0 and a strong priority con ict (SPC) if y = y0. Note that a WPC can be resolved by adding appropriate control-arcs, whereas a SPC cannot be resolved that way (see section 3.5).
Note that we have a strong priority con ict between the trivial sequence
(xk; 0) for any terminal state xk 2 Xt0erm and any sequence containing xk
as intermediate state. Such sequences are not catalytic conformal due to [
Construction of the priority con ict graph. To re ect the weak and strong priority con icts between all possible sequences resulting from X 0, we construct a priority con ict graph G = (VD [ Vterm; ED [ EW [ ES ) where the nodes correspond to sequences and the edges to priority con icts: { VD contains for all xj 2 X 0nXt0erm and the di erence vector dj = succX 0 (xi) xi, for all 2 (dj ) and all permutations of R(dj ; ) the sequence ; (xj ; dj ). { 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 re ects all strong priority con icts between sequences ; 0 stemming from distinct di erence vectors. { EW re ects all weak priority con icts between sequences ; 0 stemming from distinct di erence vectors.
The edges in ED induce a clique partition Q of G in as many cliques 3 as there are observed states in X 0 n Xt0erm resp. di erence vectors in D: VD = Q1 [ : : : [ QjDj. Moreover, each node in Vterm corresponds to a clique of size 1, so that G is partitioned into jX 0j many cliques.
Example 5. The resulting priority con ict graph G of the running example is shown in Fig. 2. 3 A clique is a subset of mutually adjacent nodes. Q1 Selection of suitable sequences. To obtain a network explaining all observations, we have to select one sequence per di erence vector dj , i.e., exactly one node from each clique Qj 2 Q. To encode the priority con icts involving terminal states, we require also to select all trivial sequences (xk; 0), i.e., all nodes from Vterm. Thus, we are interested in subsets S VD of cardinality jDj such that jS \ Qj j = 1 for each Qj 2 Q, and no SPCs occur in S Vterm.
The set of all such solutions S [ Vterm can be encoded by all vectors x 2 f0; 1gjVD[Vtermj satisfying
P 2Qj x = 1 8Qj 2 Q
x = 1 8 2 Vterm x + x 0 1 8 0 2 ES
x 2 f0; 1g 8 2 VD [ Vterm: Example 6. From G in Exp. 5, we obtain the following feasible subsets Si [ Vterm S1 = f 1(x1; d1); (x2; d2); 1(x4; d4)g, S3 = f 1(x1; d1); (x2; d2); 3(x4; d4)g, S2 = f 3(x1; d1); (x2; d2); 1(x4; d4)g, S4 = f 3(x1; d1); (x2; d2); 3(x4; d4)g. Composition of conformal networks. Every selected subset S [Vterm corresponds to a standard network PS = (P; TS ; AS ; w) which is conformal with X 0 (but not yet necessarily X 0-deterministic): { 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; { there might be weak priority con icts 0 2 EW for nodes ; 0 2 S [ Vterm which have to be resolved subsequently by inserting appropriate control-arcs. Example 7. We apply the method only to the feasible set S1 [ Vterm from Exp. 6 (all solutions for S2 [ Vterm, S3 [ Vterm and S4 [ Vterm are presented in Exp. 9). We construct the standard network presented in Fig. 3 with TS1 = fd1; d2; d4g. There is a priority con ict between (x2; d2) and (x0; 0) due to d2; 0 2 T (x2) \ T (x0). F R
PFR
Spo
R
For each of the yet reconstructed standard networks PS = (P; TS ; AS ; w) resulting from a selected subset S from the previous reconstruction step, we have to determine appropriate control-arcs in order to resolve weak priority con icts corresponding to edges 0 2 EW with ; 0 2 S [ Vterm (if any), in order to turn PS into a catalytic conformal extended Petri net PS = (P; TS ; AS ; w).
Recall that we have a weak priority con ict between two sequences and t = rt0 and intermediate states y 6= y0 with 0 if there are update vectors r 6 t; t0 2 T (y) \ T (y0) such that (y; rt) 2 but (y0; rt0 ) 2 0. This weak priority con ict has to be resolved by adding appropriate control-arcs such that { the update vector rt becomes a transition t with t 2 TA(y) but t 62 TA(y0) (or vice versa) if y; y0 62 Xterm or
0 { the update vector rt becomes a transition t with t 2 TA(y) which is disabled by control-arcs in y0 if y0 2 Xt0erm.
Inserting control-arcs This task can be done by using similar techniques as
proposed in [
Each of the so-determined control-arcs de nes a transition t with the desired properties (inheriting the standard arcs from rt and having either a read-arc or an inhibitor-arc as described above).
Remark 2. In case of a SPC involving states y = y0, the set P (y; y0) becomes empty and it is, therefore, not possible to resolve a SPC by control-arcs.
For every reconstructed standard network PS = (P; TS ; AS ; w) and any subset P 0 P containing exactly one place from P (y; y0) for every weak priority con ict, we get a catalytic conformal extended Petri net PS;P 0 = (P; TS ; AS;P 0 ; w) by inserting the respective control-arcs for all p 2 P 0. Example 8. We de ne control-arcs to resolve the WPC between (x2; d2) and (x0; 0) for the network PS1 . We obtain
P (x2; x0) = fPF R; PRg by x2 = (0; 0; 0; 1; 0)T and x0 = (0; 0; 1; 0; 0)T : Any non-empty subset of P (x2; x0) can be used to disable d2 at x0 2 Xt0erm. For PF R, x2PF R < x0PF R holds, leading to an inhibitor-arc (PF R; d2) 2 AS1;P 0 , and for PR, x2PR > x0PR holds, leading to a read-arc (PR; d2) 2 AS1;P 0 both with weight 1. The two possible alternatives are presented in Fig. 4.
F R d1
PR
d4 PFR
R
R
F R d4
PFR
d1
PR d2
Spo d2
Spo To generate the required priorities for each of the yet reconstructed extended networks PS;P 0 = (P; TS ; AS;P 0 ; w), we only need to set the priorities among all the transitions in TS according to the sequences selected for S.
Recall that every 2 S stands for a sequence =
; (xj ; dj ) = (y1; rt1 ); (y2; rt2 ); : : : ; (ym; rtm ) which induces a priority relation O indicating that the transition ti resulting from rti is supposed to have highest priority at yi. That is, 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 This implies nally that every extended network PS;P 0 = (P; TS ; AS;P 0 ; w) together with the partial order OS;P 0 constitutes an X 0-deterministic extended Petri net tting the given data X 0. Example 9. For the running example, it is left to determine the priority relations. For the extended Petri nets PS1;P 0 , we can easily verify that TAS1;P 0 (x) contains exactly one transition for all x 2 X 0, so no priorities are implied and OS1;P 0 = ; follows. For the Petri nets coming from the other sets S2; S3; S4, all possible minimal X 0-deterministic extended Petri nets are depicted in Fig. 5, 6 and 7. d2 r4
PR PFR d2 PFR
Spo d1 R d1 PR PR r3
F R F R F R F R r3 d1 PFR R d1 r4
Spo
Spo r4
R r3 r3
R r4 r3
PFR R
Spo
Spo
R r3
PFR
Spo
PFR d2 PFR
PR
PR d2
Spo To summarize, we present in this paper the steps of an integrative reconstruction method to generate all possible X 0-deterministic extended Petri nets from monotone and reproducible experimental time-series data X 0.
This approach is based on previous works for special cases: the reconstruction
of standard networks [
0
{ re ning the idea from [
Note that a preprocessing (to test the experimental data X 0 for
reproducibility and, if necessary, to handle infeasible situations) can be handled similar as in
[
In total, this integrative approach is promising for the reconstruction of networks fully tting the experimentally observed phenomena.
Our further goal is to make the new approach accessible by a suitable implementation, e.g., using Answer Set Programming as in the case of the reconstruc