The signaling petri net-based simulator: a non-parametric strategy for characterizing the dynamics of cell-specific signaling networks

Abstract

Reconstructing cellular signaling networks and understanding how they work are major endeavors in cell biology. The scale and complexity of these networks, however, render their analysis using experimental biology approaches alone very challenging. As a result, computational methods have been developed and combined with experimental biology approaches, producing powerful tools for the analysis of these networks. These computational methods mostly fall on either end of a spectrum of model parameterization. On one end is a class of structural network analysis methods; these typically use the network connectivity alone to generate hypotheses about global properties. On the other end is a class of dynamic network analysis methods; these use, in addition to the connectivity, kinetic parameters of the biochemical reactions to predict the network's dynamic behavior. These predictions provide detailed insights into the properties that determine aspects of the network's structure and behavior. However, the difficulty of obtaining numerical values of kinetic parameters is widely recognized to limit the applicability of this latter class of methods. Several researchers have observed that the connectivity of a network alone can provide significant insights into its dynamics. Motivated by this fundamental observation, we present the signaling Petri net, a non-parametric model of cellular signaling networks, and the signaling Petri net-based simulator, a Petri net execution strategy for characterizing the dynamics of signal flow through a signaling network using token distribution and sampling. The result is a very fast method, which can analyze large-scale networks, and provide insights into the trends of molecules' activity-levels in response to an external stimulus, based solely on the network's connectivity. We have implemented the signaling Petri net-based simulator in the PathwayOracle toolkit, which is publicly available at http://bioinfo.cs.rice.edu/pathwayoracle. Using this method, we studied a MAPK1,2 and AKT signaling network downstream from EGFR in two breast tumor cell lines. We analyzed, both experimentally and computationally, the activity level of several molecules in response to a targeted manipulation of TSC2 and mTOR-Raptor. The results from our method agreed with experimental results in greater than 90% of the cases considered, and in those where they did not agree, our approach provided valuable insights into discrepancies between known network connectivities and experimental observations.

Figure 1. The Model Signaling Network.

A MAPK1,2 and AKT network downstream from EGFR, which we assembled from various sources, and used for the case study analysis in this work. An edge from u to v ending with an arrow indicates an activating reaction, while an edge ending with a plunger indicates an inhibiting reaction. With the exception of TSC2, all nodes have self-inhibitory edges, which were added to model the external cellular machinery that regulates the concentration of the active form of the proteins –. Colors were selected to enhance readability of the network.

Figure 2. A High-Level Outline of the Procedure for Simulating a Signaling Network.

The input to the procedure is a signaling Petri net, S, the number of time units to simulate the network for, B, and the number of runs for which to repeat the simulation, r. The random generation of event ordering is employed to simulate the stochasticity in reaction rates and the differing times of signal arrivals.

Figure 3. The Effects of Reaction Rates on Signal Propagation.

(A) By changing the speed of signaling edge 3, the value of D at the end of a single simulation step can be reversed. If edge 3 is slower than the cascade B→C⊣D, then D will be active. If edge 3 is faster than the cascade, then D will be inactive. (B) An example of how the simulator might evaluate the individual edges during a run. In each time block, every edge is evaluated once. Each edge evaluation corresponds to one time step. Note that the order of the edge evaluation is shuffled during each time block in order to sample the space of possible relative signaling rates.

Figure 4. An Example Signaling Network and Its Corresponding Petri Net.

An example signaling network (A) and its corresponding Petri net (B). Each signaling protein in the network, A, B, and C, are designated as places p_A, p_B, and p_C. Signaling interactions become a transition node and its input and output arcs. Note that the connectivity for an activating edge differs from that of an inhibitory edge.

Figure 5. The Topological Motifs for Differing Signaling Processes.

(A) The token consumption motifs for complexing and recruitment. Transition t₁ encodes activation of v by the binding or consumption of u. Transition t₂ encodes deactivation of v by the binding or consumption of u. In both cases, the number of tokens of p_u decreases immediately after transitions t₁ and t₂ fire. (B) The token conserving motifs for PTM and GTP/ATP binding. Transition t₃ encodes enzymatic activation of v by u. Transition t₄ encodes enzymatic inhibition of v by u. In both cases, the number of tokens of p_u remains unchanged immediately after transitions t₃ and t₄ fire.

Figure 6. The Algorithm That Implements the Signaling Network Event Generator.

This routine generates the time block/firing structure. Given a set of events, E, and the number of blocks for which the SPN will be executed, n, GenerateSignalingEvents generates n blocks of events, each consisting of |E| events ordered randomly. In each block, every event in E occurs exactly once.

Figure 7. The Procedure for Simulating a Signaling Petri Net.

Simulate predicts the signal flow through the SPN S. The simulation is run for B time blocks; the results of r runs are averaged to produce the final result. Most of the work is done by the signaling Petri net execution procedure detailed in the preceding sections. This execution actually performs an individual run. This procedure takes the initial marking, m₀, and applies the sequence of transitions triggered by the event sequence, σ^e. This ordering, generated by the algorithm in Figure 6, has the dual time structure in which each block of edges contains every event in E exactly once. Each firing evaluates the effect of one transition. The markings at the end of each time block are extracted in Step 5.

Figure 8. The Algorithm for Predicting the Effect on Signal Propagation of a Targeted Manipulation.

The algorithm for predicting the effect on signal propagation of a targeted manipulation on signaling network with connectivity G. The ‘c’ and ‘p’ superscripts are used to denote parameters in the control and perturbed versions, respectively, of the SPN.

Figure 9. The Results of the TSC2 Perturbation Experiments and Simulations.

In the western blots, columns (or lanes) are as follows: (1) non-targeting (NT) control siRNA, (2) NT siRNA+EGF, (3) TSC2 siRNA, (4) TSC2 siRNA+EGF. The effect of the TSC2 siRNA on a given molecule can be assessed by comparing column 4 against column 2. For each molecule in the western blot, there is a corresponding simulation curve showing the predicted change in protein activity over time. For the purposes of this analysis, we compared the concentration change after 20 time steps (the left-most data points in the plots) for each molecule. Each simulation point corresponds to the average of 400 measurements that were computed using the procedure described in Figure 8. Experimentally derived initial states were used in the simulations. The results of both the experiments and simulations are qualitatively summarized in Table 3.

Figure 10. The Predicted Response of the Network to an mTOR-Raptor Perturbation.

The predicted response of the network to a mTOR-Raptor perturbation in the (A) MDA231 and (B) BT549 cell lines. Our method predicts that the amount of available AKT increases in response to the perturbation, which is in agreement with results published in the literature ,. Our method also predicts that the activity-level of p70S6K in the MDA231 cell line decreases in response to the perturbation, which has been observed experimentally . Each point corresponds to the average of 400 measurements that were computed using the procedure described in Figure 8.