Multilevel model for the investigation of oncoantigendriven vaccination effect
 Francesca Cordero^{1}Email author,
 Marco Beccuti^{1},
 Chiara Fornari^{1, 2},
 Stefania Lanzardo^{2},
 Laura Conti^{2},
 Federica Cavallo^{2},
 Gianfranco Balbo^{1, 3} and
 Raffaele Calogero^{2}
https://doi.org/10.1186/1471210514S6S11
© Cordero et al.; licensee BioMed Central Ltd. 2013
Published: 17 April 2013
Abstract
Background
Cancer stem cell theory suggests that cancers are derived by a population of cells named Cancer Stem Cells (CSCs) that are involved in the growth and in the progression of tumors, and lead to a hierarchical structure characterized by differentiated cell population. This cell heterogeneity affects the choice of cancer therapies, since many current cancer treatments have limited or no impact at all on CSC population, while they reveal a positive effect on the differentiated cell populations.
Results
In this paper we investigated the effect of vaccination on a cancer hierarchical structure through a multilevel model representing both population and molecular aspects. The population level is modeled by a system of Ordinary Differential Equations (ODEs) describing the cancer population's dynamics. The molecular level is modeled using the Petri Net (PN) formalism to detail part of the proliferation pathway. Moreover, we propose a new methodology which exploits the temporal behavior derived from the molecular level to parameterize the ODE system modeling populations. Using this multilevel model we studied the ErbB2driven vaccination effect in breast cancer.
Conclusions
We propose a multilevel model that describes the interdependencies between population and genetic levels, and that can be efficiently used to estimate the efficacy of drug and vaccine therapies in cancer models, given the availability of molecular data on the cancer driving force.
Keywords
Background
Systems biology is increasingly used to get insights into the functioning of complex biological networks. Specifically, the use of mathematical formalisms to investigate the mechanisms affecting tumor growth and maintenance upon vaccination or drug treatment might represent a powerful instrument to efficiently guide the design of long and expensive in vivo experiments [1].
Building network models that accurately represent either biochemical pathways, celltocell interactions, or regulation networks is necessary for different purposes. Indeed, a model provides the basis for a clear description of the interactions involved in a biological system. However, to be useful, a model must be precise and suitable for an analysis that helps in getting a better understanding of the phenomenon under investigation and appropriate formalisms must be used to achieve this goal. Furthermore, when the objective of the study is the behavior of a biological system described at the level of a biochemical reaction scheme, the completion of the modelling process sets the ground for a sensitivity analysis of the model where, at the level of molecule concentrations, it is possible to perturb the net representation or the reaction rates to study the influence of specific elements of the network on the overall functionality of the system. From a structural point of view, a qualitative analysis of the model can be used to select key elements that may suggest interesting features of the experimental system which are worth of detailed investigations (e.g. therapeutic targets). All these points highlight the need of a strong integration between computational modeling and quantitative experimental data. The paper by Kreeger and Lauffenburger [2] reports examples of recently proposed integrations of this type in the field of cancer systems biology. However, even thought a pathwaycentric approach is widely and successfully used to investigate cancer in terms of molecular effects, when a specific gene or protein is identified to make a contribution to pathology, it is not easy to determine how its influence is propagated at population level, unless the interaction between molecular effects and population dynamics is specifically addressed by the model.
In this paper we propose a new approach, which allows to describe in a single multilevel model different dynamics levels (e.g. molecular regulation network and cell population) of a complex biological system providing a way of highlighting the interactions among different levels and making easier the model parameter definition. For instance the difficulty of defining the parameters of a population model, in the presence of few measured results, can be overcome deriving such parameters directly by a molecular network that mimics the most relevant biochemical reactions occurring into a cell population, thus accounting for the presence of environmental changes, mutation, and noise in intracellular biochemical reactions.
Even if it is true that with this methodology, dealing with a lack of information is only moved one level down in the modeling process, (functionally) simpler parameters are now used in the basic model of our case study, leaving to the solution of the molecular model the burden of deriving the complex representations of the proliferation parameters that are now allowed to be expressed with intricate functions of time.
As we just said, our case study has been modeled with a 2level representation. The first level describes regulation aspects of proliferation considering gene interactions; specifically, our experimental model of carcinogenesis is driven by ErbB 2 [3]. This molecular network is designed using the Petri Net [4] formalism which is quite suitable to build models of this type and which allows to compute qualitative and quantitative properties of the experimental system with numerical and analytical methods. Moreover, PNs offer the possibility of representing a reaction scheme as a graphical diagram that supports the comprehension of the behavior of the real system with simple to understand, yet precise descriptions.
The second level describes the population interactions in the ErbB 2driven carcinogenesis, and is based on the model presented in Fornari's paper [5], where a system of ODEs was used to describe the progression of malignant tumors, assuming the validity the CSC theory [6]. Our ODE model takes into account the main properties of CSCs: tumorigenic capacity, selfrenewal, and differentiation into nonstem cells. The hierarchical organization of the tumor is guaranteed both from the growth and progression as well as from the differentiation capacity characterizing CSC subpopulations. CSCs give rise to committed Progenitor Cells (PCs) characterized by a rapid proliferation rate. PCs are able to completely differentiate into Terminally differentiated Cells (TCs). In the paper by Fornari et al., parameters characterizing the behavior of proliferation, death, and differentiation of tumor cell populations are assumed to be affected by external events such as vaccination or pharmacological treatments and are tuned using experimental data based on the tumor mass growth trend observed in mice after a subcutaneous injection of cancer cells. However, to make such cellular model interesting for the biological community, it is necessary to link proliferation and differentiation parameters to the molecular events which control them, thus allowing the visualization at cellular level of perturbations made on the underlying molecular network.
To investigate how the perturbation at molecular level impacts on population models, we used a simulative approach to show the effects of well known inhibition of progression of multifocal preneoplastic lesions [7] in ErbB2driven carcinogenesis, by means of chronic vaccination and we provide an example of new hypothesis that can be generated using such models and subsequently validated with biological experiments.
Results and discussion
In this section, we first discuss the new proposed approach in details, and then we show how it can be used to study and analyze the effects of vaccination on a carcinogenesis driven by ErbB 2 [3] receptor family considering both population and molecular aspects.
Workflow
Model definition
In this first step the biological system is represented by a multilevel model, where the number of levels is chosen according to the phenomenon under study. As already highlighted in the background Section, focussing at each level on different aspects of the problem under study, model creation and parameterization are made easier. For instance, our case study which is concerned with the carcinogenesis driven by ErbB 2, was modeled by 2level model, where the former level describes the molecular regulatory network and the latter one the cell populations. In details, the dynamics of the low level, namely at the gene and molecular scale, is modeled using the PN formalism. It provides an explicit and intuitive representation of the signaling cascade controlled by the ErbB receptors family, capturing the relevant biochemical reactions involved in the regulation aspects of proliferation. Instead, the dynamic of the high level system, namely at the cell population scale, which describes the interactions between different subpopulations of cells is represented by a system of ODEs, following a trend that has already been established in the literature [8].
Model consistency and correctness
This second step that naturally follows the model definition phase is focused on the validation and verification of the accuracy and correctness of the representation. Among several methodologies, used to perform these analysis, we exploit structural validation and model checking. The model structure is validated using the Psemiflow analysis, and thus identifying the set of places where a given kind of correlated matter is preserved during the evolution of the model in order to make sure that mass conservation laws are respected. On the other hand, model checking is used to verify the consistency and correctness of the model with respect to wellknown properties found in the literature and expressed through Computational Tree Logic (CTL), a type of temporal logic. In our case study, model checking is used to verify if the growth factors stimulation always leads to a production of specific protein complexes. Details on how model checking is used in our case study also to analyze other properties of the model will be provided at some length in the Methods Section of this paper.
Multi level model interactions
After the creation and the validation of models, it is necessary to define how models interact. In our case study, to link populations proliferation parameters with regulation events we identify for each level a set of interaction points. Specifically, for the molecular level we select a set of (PN) places playing a pivotal role in cell proliferation corresponding to Bad, cyclinD, and NFkB proteins. On the other hand, at population level, we select proliferation rates of CSC and PC as the parameters that mostly depend on biochemical reaction dynamics. Hence, we specify interactions defining each proliferation rate as the product of three functions representing the temporal behaviors of protein targets.
Model dynamics
The last step is related to the analysis of the global model dynamics. First, the corresponding ODE system is automatically derived from the PN model.
In general, this system of ODEs can be very large and complex, thus a preliminary reduction phase is performed to obtain a suitable system of ODEs. Specifically this phase consists in a downsizing of the ODE number by identifying those equations which are redundant using PN structural properties such as Psemiflows. Indeed, we derive a set of ODEs from the PN model and, then, for each minimal Psemiflow of the system one equation from the ODE system is removed. After having reduced the complexity of the model, the temporal dynamics of the quantities contained in the places which play a pivotal role in cell proliferation are studied through numerical integration of the derived ODE system. Hence, the obtained quantities are used as parameters in the ODE system modeling the cell proliferation, and it is solved by numerical integration.
Referring to our case study, in the following paragraphs we show how our methodology can be be put in practice discussing first the individual components of our twolevel model.
First level: molecular regulation network
The main effects of the downstream signaling is the production of Pip_{3} that leads to the activation of Akt, as reported in portion B of the network. The production of Pip_{3}, which is a second messenger involved in the regulation of different processes, is catalyzed by PI3K starting from Phosphatidylinositol 4,5triphosphate, Pip_{2}. In portion B a set of reactions involved in the regeneration of Pip_{2} is also reported. Its recovery results from the contribution of the Ptendependent dephosphorylation of Pip_{3}.
With respect to Birtwistle's work, we extend the network with three additional blocks (Portions C, D, and E). Portion C describes the downstream effects of Akt activation. Akt has a critical regulatory role in many cellular processes, and in particular in cancer progression. As described before, we decided to focus the effects of Akt on three targets:

♦ the transcription factor Bad the proliferation action of Akt is mediated through the direct inhibition of this proapoptotic signal,

♦ the activation of cyclinD Akt occurs at the G1S transition of the cell cycle via phosphorylation and inhibition of glycogen synthase kinase 3beta (GSK3β) that stabilizes cyclin D1,

♦ the transcription factor nuclear factorkappa B, NFkB Akt promotes NFkB activity since it directly phosphorylates IkappaB kinase α, IKKα, to activate NFkB whose broad oncogenesis activity  through its ability to control cell proliferation and to suppress apoptosis  is well known.
Another important regulation of cell growth by Akt regards its primary effect on mTOR whose action is depicted in portion D of the network. mTOR is associated with two complexes: the rapamycinsensitive TORC1 complex which controls S6K phosphorylation and 4EBP1 to regulate translation, and TORC2 that controls the phosphorylation of Akt. The activation of TORC1 by Akt involves the phosphorylation of TSC2, which reveals a negative regulatory effect on mTOR controlled by the GTPase Rheb.
Finally, portion E specifies the cascade of TLR2. Functional analysis of mammalian TLRs has revealed that they recognize specific patterns involved in the cell proliferation. The signaling pathway via TLR2 recruits the adapter protein MyD88. Upon stimulation, MyD88 recruits IL1 receptorassociated kinase (IRAK) to TLR2. IRAK is activated by phosphorylation and then associated with TRAF6, leading to the activation of two distinct signaling pathways, and finally to the activation of JNK and NFkB.
Overall, this network is a modification of that proposed by Birtwistle in order to account the characteristics of preclinical breast cancer model based on BALB/c mice transgenic for the transforming rat ErbB2 oncogene, BALBneuT mice. BALBneuT mice develop breast cancer with 100% penetrance [10]. These animals are transgenic for a mutated ErbB2 rat gene, encompassing a single point mutation that replaces the valine residue at position 664 in the transmembrane domain with glutamic acid favoring ErbB2 homodimerization thus transforming the ErbB2 protooncogene into a dominant transforming oncogene [11]. In vivo experiments have shown that PI3K represents an important element in the ErbB2 signal transduction since antiErbB2 antibodies impair PI3K/Aktmediated tumorigenic effects [12], these experiments also demonstrate the ability of ErbB 2 to activate directly Akt without the involvement of growth factors. Moreover, the choice of adding the TLR2 contribution to the proliferation pathway derives from the observation that the TLR2 receptor shares the PI3K activation network [13] with ErbB2 [14], and accounts for recent results that show TLR2 to be expressed by breast cancer cell lines [15] and to be involved in cancer invasiveness.
Second level: cell population model
where N_{ CSC }, ${N}_{P{C}_{i}}$ and N_{ TC } are the numbers of cancer stem cells, progenitor cells, and terminal cells, respectively. Notice that the terms characterizing these equations depend on 4 parameters: ω_{ CSC } and ω_{ PC } describe the proliferation rates; γ_{ PC } represents the bidirectional interconvertibility parameter that involves CSC, PC_{1}, PC_{2} and PC_{3} subpopulations; d_{ i } indicates the death rate; and η_{ i } describes the differentiation rates.
Parameter definition
TLR2 and Neu expression on TuBo (E) and serial passages mammospheres (P1, P2, P3).
Cell line  ErbB2 (%)  TLR2 (%) 

E  52.7 ± 4.2  8 ± 1 
P1  51.9 ± 12.7  21 ± 3 
P2  45.5 ± 10.5  24.5 ± 5.3 
P3  41.2 ± 10.4  27.1 ± 3.8 
The initial marking of three growth factors is defined as a function which models injections happening at regular time intervals. For what concern the parameters at population level, the TCs/CSCs ratio was defined indirectly on the basis of spheroids that can be produced starting from a cell culture of 1000 TUBO cells. Since spheroids are clonal, i.e. each spheroid is derived from a unique CSC, counting the number of spheroids allows to quantify the number of CSCs present in the starting cell culture (15 +/5 CSCs every 1000 TUBO cells). The differentiation, death, and bidirectional interconvertibility parameters are set as reported in Fornari's work [5]. Otherwise, proliferation rates of cell population are defined considering several proteins dynamics at regulation level.
Model consistency and correctness
The consistency and correctness of the model has been verified applying three preliminary checks. The first check is based on Psemiflows that, as explained in details in the Methods Section, can be used to identify the sets of places where a given kind of correlated matter is preserved. In this way using the biological reactions involved in the model construction it is possible to identify sets of places that must appear in the same Psemiflows.
we observe that the places representing the conservation of Pten are correctly in the same Psemiflow: Pten, Pip_{3}:Pten, Pip_{2}:Pten, Pten:Pip_{2}:Pip_{3}.
The second check is based on model checking: a technique that provides a useful quality control for the development of large scale complicated models. As described in the Methods Section, model checking is based on the use of temporal logic to specify system behavioral proprieties which can be processed automatically using computationally efficient procedures that determine whether they are effectively reproduced by a model. In this work we have used model checking to verify the consistency and correctness of our model with respect to a set biological reactions used for its construction, and to other wellknown properties found in literature.
Finally the last check consists in verifying that the quantitative behavior of the model is consistent with the literature, for instance according to results presented in [18] we verified that Pten inhibition leads to an inhibition of the proliferation pathway.
Multilevel model interaction
To link proliferation (population level) with regulation events (regulatory network) it is necessary to set interaction points from both models. Those points of interactions are selected through the singling out of proteins (at regulatory level) directly involved in the phenomenon under study (at population level). As hinted before, for what concerns the population level, we select proliferation parameters, i.e. ω_{ CSC } and ω_{ PC } . Instead, at the molecular level we select proteins which have a pivotal role in cell proliferation, i.e. cyclin D, NFkB and BAD. The interaction is then defined assigning at proliferation parameters specific values deduced from those target proteins. In detail, three functions representing the temporal behaviors of cyclinD, NFkB and BAD are created for both CSC and PC regulatory networks. These functions are obtained from the solution of the ODE systems corresponding to the (first level) molecular network of CSC and PC. Proliferation rates are then evaluated as the product of the three functions, which take different values in CSC and in PC regulatory networks.
The vaccination backlash, applied at molecular level, is directly reflected on protein targets.
Model dynamics
In this last step we describe the two experimental analyses performed within our case study.
Effects of ErbB2 vaccination
Involvement of TLR2 in CSC proliferation
We are presently evaluating with in vitro experiments, the effect of TLR2 silencing in TUBO cell proliferation to confirm these simulation results. It is notable that literature data [16] indicate that TLR2 mediate innovation by the activation of the NFkB pathway. This finding together with our observation that TLR2 positive cells are mainly associated with subpopulation of cancer cells enriched for CSC (Table 1), suggest that TLR2 might play some important role in CSC invasiveness. Thus, the TLR2 network might represent an interesting starting point to design a network controlling the parameters linked to CSC from PC differentiation, since invasiveness is associated with undifferentiated cells, i.e. CSC, and is lost in fully differentiated cells, i.e. TC.
Conclusions and perspectives
In this paper we propose a novel approach in which a multilevel model is constructed and where molecular networks are used to estimate certain parameters of a cell population model based on a system of ODEs. With this method we have been able to reproduce at a qualitative level the effect of antiErbB2 chronic vaccination in BALBneuT model. Although the model needs some refinement to provide a punctual representation of vaccination, i.e. aligning the time line of the computational model with in vivo data, it successfully supports the idea that new in vitro/in vivo experiments can be designed to test hypothesis that are formulated on the basis of the solution of the model. Furthermore, our approach can be extended to consider the immunological tumor microenvironment by adding new equations in the ODE system of the population representation and by defining their parameters on the basis of a celltocell network, instead of a genetic network. This might be particularly interesting in the area of combined treatment development. Tumor vaccination alone is not sufficient to eradicate the disease, but combined with other immunopharmacological treatments, affecting the CSC differentiation rate might represent an interesting approach in the area of tertiary cancer prevention, i.e. reducing the negative impact of disease by restoring functions and reducing diseaserelated complications.
Methods
The following section reports the details of the biological techniques used for the experiments as well as the notation and the basic definitions of formalism and algorithms used for the analysis discussed in this paper. Part of this section is similar to an analogous description reported in our previous work [23] and is included here only to help the reader in understanding the subsequent discussion.
Flow cytometry analysis
E, P1, P2, and P3 cells were collected after 7 d of culture and disaggregated through enzymatic and mechanical dissociation. They were then washed in PBS (SigmaAldrich) supplemented with 0.2% BSA and 0.01% sodium azide (SigmaAldrich) and stained for the membrane antigen Toll Like Receptor (TLR) 2 using an Alexa Fluor647conjugated antiTLR2 and for ErbB2 antigen using Ab4 moAb (Calbiochem) followed by FITC antimouse IgG (Dako Cytomation) as secondary Abs. Samples were collected and analyzed with a CyAn ADP Flow Cytometer and Summit 4.3 software (DakoCytomation).
Petri Net formalism
A brief introduction of the PN formalism [24] is provided here to help the reader in understanding the model of the proliferation pathway that we have used in this work.
PNs have been first proposed for the representation of biological pathways by Reddy et al. in [4]. Subsequently, many other researchers discussed the advantages of using PNs to model biological systems [25] because of their ability of representing reaction systems in a natural graphical manner and of their capability of allowing the computation of qualitative and quantitative information about the behavior of these systems. Most of the analysis of PN models in Systems Biology is performed using simulative approaches [26, 27] where transitions occur after random firing delays. The extension of the basic PN formalism with stochastic firing delays of this type has been proposed in the literature with the definition of the socalled Stochastic Petri Nets (SPNs) [28–30] which allow the automatic construction of an underlying Continuous Time Markov Chain (CTMC) that can be studied using numerical or simulative techniques and that can also be translated into systems of ODEs when average results are sufficient for the analysis [31]. Several papers recently published by Heiner et al. [32, 33] show that the available biological data can be analyzed by means of PNs formalism in order to obtain the general behavioral patterns and timing measures of biological entities. In general, biological models are affected by a high level of complexity due to the functional properties of the systems that are considered. The interaction of qualitative and quantitative analysis is necessary to check a model for consistency and correctness as we will show in the rest of this paper using PNs. In details, PNs are bipartite directed graphs with two types of nodes: places and transitions. The places, graphically represented as circles, correspond to the state variables of the system (e.g. enzymes, compounds, etc.), while the transitions, graphically represented as boxes, correspond to the events (e.g. interactions among biochemical entities) that can induce state changes. The arcs connecting places to transitions (and vice versa) express the relations between states and event occurrences. Places can contain tokens (e.g molecules of the corresponding entities) drawn as black dots within places. The state of a PN, called marking, is defined by the number of tokens in each place. The evolution of the system is given by the firing of enabled transitions, where a transition is enabled if only if each input place contains a number of tokens greater or equal than a given threshold defined by the cardinality of the corresponding input arc. A transition occurrence/firing removes a fixed number of tokens from its input places and adds a fixed number of tokens into its output places (according to the cardinality of its input/output arcs).
The set of all the markings that the net can reach, starting from the initial marking through transition firings, is called the Reachability Set (RS). Instead, the dynamic behavior of the net is described by means of the Reachability Graph (RG), an oriented graph whose nodes are the markings in the RS and the arcs represent the transition firings that produce the corresponding marking changes.
Here we recall briefly the notation and the basic definitions that are used in the rest of the paper.
Definition: Petri Net
A PN system is a tuple $\mathcal{N}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{0.3em}{0ex}}\left(P,T,{I}^{},{I}^{+},{\mathsf{\text{m}}}_{0}\right)$ where:

♦ P = {p_{ i }} is a finite and non empty set of places.

♦ T = {t_{ i }} is a finite and non empty set of transitions with P ∩ T = ∅.

♦ ${I}^{},{I}^{+}:T\times P\to \mathbb{N}$are the input, output, that define the arcs of the net and that specify their multiplicities.

♦ ${\mathsf{\text{m}}}_{0}:P\to \mathbb{N}$is a multiset on P representing the initial marking of the net.
A marking m (or state) of a PN is a multiset on P. A transition t is enabled in marking m iff I^{  }(t, p) ≤ m(p), ∀p ∈ P , where m(p) represents the number of tokens in place p in marking m. Enabled transitions may fire, so that the firing of transition t in marking m yields a marking m' = m + I^{+}(t)  I^{}(t). Marking m' is said to be reachable from m because of the firing of t and is denoted by m[t〉m'. The firing of a sequence σ of transitions enabled at m and yielding m' is denoted similarly: m[σ〉m'.
How to check model consistency and correctness
Before describing how it is possible to study the temporal dynamics of a phenomenon described with a PN model, we discuss two preliminary analysis steps useful to verify the model consistency and correctness, based on Psemiflows and CTL properties. The first step can be used to identify the set of places where a given kind of correlated matter is preserved during the evolution of the model, while the latter allows the modeler to verify more complex behavioral properties.
Mathematically, Psemiflow [34] can be defined as follows:
Definition: Psemiflow
Given a Petri Net, let C be the Incidence Matrix whose generic element ${c}_{t,p}={I}^{+}\left(t,p\right){I}^{}\left(t,p\right)$describes the effect of the firing of transition t on the number of tokens in the place p; and let $x\in {\mathcal{Z}}^{\leftP\right}$be a place vector; then a Psemiflow is a place vector x such that it represents an integer and nonnegative solution of the matrix equation xC = 0.
All the Psemiflows of a PN can be expressed as linear combinations of a set of minimal Psemiflows, and the support of a Psemiflow $\mathcal{F}\phantom{\rule{0.1em}{0ex}}$, denoted supp$\left(\mathcal{F}\right)$ can be defined as the set of nodes corresponding to the nonzero entries of $\mathcal{F}\phantom{\rule{0.1em}{0ex}}$. Using supp$\left(\mathcal{F}\right)$, each Psemiflow $\mathcal{F}\phantom{\rule{0.1em}{0ex}}$ allows the computation of a corresponding weighted sum of tokens contained in a subset of places of the net that remains constant through the entire evolution of the model; this constant ia called Pinvariant.
In a biological context, where tokens represent compounds, enzymes etc., the interpretation of such Pinvariant is relatively simple: the places of supp$\left(\mathcal{F}\right)$ represent the portion of the PN where a given kind of correlated matter is preserved. Obviously when all the places of a net belong to at least one Psemiflow, then the markings of the places are bounded and the state space of the net is finite.
Finally it is important to observe that Psemiflow analysis involves only the structural proprieties of the net and is thus independent of the initial marking of the PN.
The second preliminary analysis step allows the user to verify more complex behavioral properties using the model checking technique, a wellestablished formal method that is widelyused for ascertaining the correctness of reallife systems. It requires a description of the system usually given in some high level modeling formalism such as PN, and the specification of one or more desired properties of the system, normally using temporal logics. Given this input, the model checker can derive the system behavior (e.g. generating the system RG) and automatically verify whether or not each property is satisfied though a systematic and exhaustive exploration of the RG.
Among the many available temporal logic formalisms, we choose the Computational Tree Logic, CTL, a branchingtime logic that extends propositional logic used for describing states, with operators for reasoning over time and nondeterminism.
In details the following temporal operators are considered in CTL: Xp meaning that the proposition p is true at the next transition, Gp meaning that p is always true, Fp meaning finally true, pUq meaning that p is true until q becomes true. For reasoning about nondeterminism, the two following path quantifiers are used: Ap meaning that p is true on all paths and Ep meaning that p is true on some path. All the temporal operators have to be immediately preceded by a path quantifier, hence AXGp is not a valid CTL formula, since the temporal operator G is not preceded by a path operator. Moreover, atomic propositions consist of statements on the current token situation in a given place, and they can be recursively composed into more complex propositions using the standard logical operators: ¬ (negation), ∧ (conjunction), ∨ (disjunction), and ⇒ (implication). Hence the CTL can be sufficiently expressive to encode a wide range of biological queries:

♦ reachability queries: there is a cascade of reactions that lead to the production of a protein p  EFp;

♦ pathway queries: an enzyme can reach an activation state s through a substratebound state sb  EF (sb ∧ EF (s)); a cell can reach a state s without violating a certain constrain c  E(cUs); a protein p can be synthesized without a set of transcriptional factors q  E(¬qUp);

♦ steady state query a certain state s of a network is a steady state  s ⇒ EG(s); an enzyme can stay in active or inactive state  EF (AGs); an enzyme exhibits a cyclic behavior with respect to the presence of an activator or inhibitor z  EG((z ⇒ EF ¬z) ∧ ((¬z ⇒ EFz).
Finally, it is worth noting that, while CTL is an extremely powerful and flexible language to describe special properties, it can be used only by skilled users since it requires a certain experience to correctly express the specification of the desired behavior. Furthermore model checking technique, in case of very complex systems, can require substantial computational resources (in terms of memory and time) since it needs the generation of the RG of the system. In literature several approaches based on efficient RG encoding and manipulation [35] or level concept and monotonic liveness [33] have been proposed to cope with such problem.
How to analyze temporal dynamics of the modeled system
To model and study the temporal dynamics of a PN we have to introduce temporal specification in the formalism. As we already said at the beginning of this Section, the most common timed extension of PN is Stochastic PN, SPN [28] in which exponentially distributed random delays (interpreted as durations of certain activities) are associated with the firings of the transitions. In details a SPN can be defined as a pair $\left(\mathcal{N},w\right)$, where $\mathcal{N}\phantom{\rule{0.1em}{0ex}}$ is a Petri net and $w:{\mathbb{N}}^{P}\times T\to \mathbb{R}$ is a (possibly marking dependent) function that assigns to each transition of the net the rate of a negative exponential distribution of the firing delay.
Hence, for any transition t it is necessary to specify a function w(m, t), so that when t is enabled in a marking m then w(m, t) has to be evaluated to provide the rate of t in m.
where E is the number of interacting entities and X_{ i }(t) represents the amount of the i^{ th } entity at time t, N_{ i } the number of reactions in which the i^{ th } entity is involved, the parameters w_{ ij } the rate describing the speeds of these reactions, and the parameters g_{ ijh } the socalled kinetics orders which depend on the stoichiometry and on the mechanisms of the reactions. The ODEs and the initial amount of the different entities can be automatically obtained from the SPN representation, and numerical integration of the ODEs is performed to calculate the quantities at a given time instant. It is important to observe that when the number of tokens increases the quantitative behavior obtained applying the stochastic approach tends to that obtained from the ODEs [31]. Hence from the SPN description of the biological phenomenon, the choice of using one of the two approaches (stochastic or deterministic) for studying the behavior of the system is left to the analyst who decides on the basis of the objectives of his/her study. In this paper we use the deterministic approach because it allows a faster and simpler evaluation of the proposed pathway.
How to use the Psemiflow to reduce the number of ODEs
Psemiflows can be used to reduce the number of ODEs representing the behavior of the system, by identifying those which are redundant. Indeed, as already explained, Psemiflows, can be used to derive the set of places in which the total mass is preserved so that the sum of their corresponding ODEs yields a zero identity. Hence, for each minimal Psemiflow in the model we can select one place belonging to it and rewrite its corresponding variable as a linear combination of the other variables in the same Psemiflow. In this way we reduce the number of ODEs of one unit for each minimal Psemiflow present in the system.
Declarations
Acknowledgements
This study was funded under the auspices of European Consortium for Anticancer Antibody Development (EUCAAD) 200755, and by grants from the Italian Association for Cancer Research (AIRC IG 11675), the Epigenomics Flagship Project EPIGEN, MIURCNR, the Italian Ministero dell'Università e della Ricerca, and the University of Torino and Regione Piemonte. We appreciate the support by King Abdulaziz University of Saudi Arabia. We thank the reviewers for their insightful suggestions.
Declarations
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 6, 2013: Selected articles from the 10th International Conference on Artificial Immune Systems (ICARIS). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S6.
Authors’ Affiliations
References
 Hlavacek W: How to deal with large models?. Molecular Systems Biology. 2009, 5:Google Scholar
 Kreeger P, Lauffenburger D: Cancer systems biology: a network modeling perspective. Carcinogenesis. 2010, 31: 28. 10.1093/carcin/bgp261.PubMed CentralView ArticlePubMedGoogle Scholar
 Pannellini T, Forni G, Musiani P: Immunobiology of her2/neu transgenic mice. Breast Disease. 2004, 20: 3342.PubMedGoogle Scholar
 Reddy V, Mavrovouniotis M, Liebman M: Petri net representation in metabolic pathways. Proc Int Conf Intelligent Systems for Molecular Biology. 1993, 328336.Google Scholar
 Fornari C, Cordero F, Manini D, Calogero R, Balbo G: Mathematical approach to predict the drug effects on cancer stem cell models. Proceedings of the CS2Bio 2nd International Workshop on Interactions between Computer Science and Biology. Reykjavik, Iceland. 2011Google Scholar
 Hemmings C: The elaboration of a critical framework for understanding cancer: the cancer stem cell hypothesis. Pathology. 2010, 42: 105112. 10.3109/00313020903488773.View ArticlePubMedGoogle Scholar
 Quaglino E, Iezzi M, Mastini C, Amici A, Pericle F, Carlo ED, Pupa S, Giovanni CD, Spadaro M, Curcio C, Lollini P, Musiani P, Forni G, Cavallo F: Electroporated DNA vaccine clears away multifocal mammary carcinomas in her2/neu transgenic mice. Cancer Research. 2004, 64: 28582864. 10.1158/00085472.CAN032962.View ArticlePubMedGoogle Scholar
 Birtwistle M, Hatakeyama M, Yumoto N, Ogunnaike B, Hoek J, Kholodenko B: Liganddependent responses of the ErbB signaling network: experimental and modeling analyses. Molecular Systems Biology. 2007, 3: 144160.PubMed CentralView ArticlePubMedGoogle Scholar
 Liang J, Slingerland J: Multiple roles of the PI3K/PKB (Akt) pathway in cell cycle progression. Cell Cycle. 2003, 2: 339345.View ArticlePubMedGoogle Scholar
 Boggio K, Nicoletti G, Carlo ED, Cavallo F, Landuzzi L, Melani C, Giovarelli M, Rossi I, Nanni P, Giovanni CD, Bouchard P, Wolf S, Modesti A, Musiani P, Lollini P, Colombo M, Forni G: Interleukin 12mediated prevention of spontaneous mammary adenocarcinomas in two lines of Her2/neu transgenic mice. J Exp Med. 1998, 188: 589596. 10.1084/jem.188.3.589.PubMed CentralView ArticlePubMedGoogle Scholar
 Gullick W, Bottomley A, Lofts F, Doak D, Mulvey D, Newman R, Crumpton M, Sternberg M, Campbell I: Three dimensional structure of the transmembrane region of the protooncogenic and oncogenic forms of the neu protein. EMBO journal. 1992, 11: 4348.PubMed CentralPubMedGoogle Scholar
 Porzia A, Lanzardo S, Citti A, Cavallo F, Forni G, Santoni A, Galandrini R, Paolini R: Attenuation of PI3K/Aktmediated tumorigenic signals through PTEN activation by DNA vaccineinduced antiErbB2 antibodies. Journal of Immunology. 2010, 15: 41704177.View ArticleGoogle Scholar
 Li L: Regulation of innate immunity signaling and its connection with human diseases. Current Drug Targets Inflamm Allergy. 2004, 3: 8186. 10.2174/1568010043483863.View ArticleGoogle Scholar
 Ghatak S, Misra S, Toole B: Hyaluronan constitutively regulates ErbB2 phosphorylation and signaling complex formation in carcinoma cells. Journal Biology Chemical. 2005, 280: 88758883.View ArticleGoogle Scholar
 Xie W, Huang Y, Xie W, Guo A, Wu W: Bacteria peptidoglycan promoted breast cancer cell invasiveness and adhesiveness by targeting tolllike receptor 2 in the cancer cells. PLoS One. 2010, 26: 115.Google Scholar
 Rovero S, Amici A, Carlo ED, Bei R, Nanni P, Quaglino E, Porcedda P, Boggio K, Smorlesi A, Lollini P, Landuzzi L, Colombo M, Giovarelli M, Musiani P, Forni G: DNA vaccination against rat her2/Neu p185 more effectively inhibits carcinogenesis than transplantable carcinomas in transgenic BALB/c mice. Journal of Immunology. 2000, 165: 51335142.View ArticleGoogle Scholar
 O'Reilly K, Rojo F, She Q, Solit D, Mills G, Smith D, Lane H, Hofmann F, Hicklin D, Ludwig D, Baselga J, Rosen N: mTOR inhibition induces upstream receptor tyrosine kinase signaling and activates Akt. Cancer Research. 2006, 66: 15001508. 10.1158/00085472.CAN052925.PubMed CentralView ArticlePubMedGoogle Scholar
 Hers I, Vincent E, JM JT: Akt signalling in health and disease. Cell Signalling. 2011, 23: 15151527. 10.1016/j.cellsig.2011.05.004.View ArticlePubMedGoogle Scholar
 Drebin J, Link V, Stern D, Weinberg R, Greene M: Downmodulation of an oncogene protein product and reversion of the transformed phenotype by monoclonal antibodies. Cell. 1985, 41: 697706.View ArticlePubMedGoogle Scholar
 Katsumata M, Okudaira T, Samanta A, Clark D, Drebin J, Jolicoeur P, Greene M: Prevention of breast tumour development in vivo by downregulation of the p185neu receptor. Nature Methods. 1995, 1: 644648. 10.1038/nm0795644.Google Scholar
 Klapper L, Vaisman N, Hurwitz E, PinkasKramarski R, Yarden Y, Sela M: A subclass of tumorinhibitory monoclonal antibodies to ErbB2/HER2 blocks crosstalk with growth factor receptors. Oncogene. 1997, 14: 20992109. 10.1038/sj.onc.1201029.View ArticlePubMedGoogle Scholar
 Xu F, Lupu R, Rodriguez G, Whitaker R, Boente M, Berchuck A, Yu Y, DeSombre K, Boyer C, Bast R: Antibodyinduced growth inhibition is mediated through immunochemically and functionally distinct epitopes on the extracellular domain of the cerbB2 (HER2/neu) gene product p185. Int J Cancer. 1993, 53: 401408. 10.1002/ijc.2910530310.View ArticlePubMedGoogle Scholar
 Cordero F, Horvath A, Manini D, Napione L, Pierro MD, Pavan S, Picco A, Veglio A, Sereno M, Bussolino F, Balbo G: Simplification of a complex signal transduction model by the application of invariants and flow equivalent server. Theoretical Computer Science. 2011, 412: 60366057. 10.1016/j.tcs.2011.06.013.View ArticleGoogle Scholar
 Peterson J: Petri Net Theory and the Modeling of Systems. 1981, Upper Saddle River, NJ, USA: Prentice Hall PTRGoogle Scholar
 Hofestädt R: A Petri Net Application of Metabolic Processes. Journal of System Analysis, Modeling and Simulation. 1994, 16: 113122.Google Scholar
 Goss P, Pecoud J: Quantitative modeling of stochastic systems in molecular biology by using stochastic Petri nets. Proc Natl Acad Sci. 1998, 95 (12): 67506755. 10.1073/pnas.95.12.6750.PubMed CentralView ArticlePubMedGoogle Scholar
 Hofestädt R, Thelen S: Quantitative modeling of biochemical networks. In Silico Biology. 1998, 1 (6):Google Scholar
 Molloy MK: Performance Analysis using Stochastic Petri Nets. IEEE Transaction on Computers. 1982, 31 (9): 913917.View ArticleGoogle Scholar
 Natkin S: Les réseaux de Petri stochastiques et leur application à l'évaluation des systèmes informatiques. Thèse de Docteur Ingégneur, CNAM. 1980Google Scholar
 Balbo G: Introduction to Stochastic Petri Nets. Lectures on Formal Mathods and Performance Analysis, Volume 2090 of LNCS. Edited by: Brinksma E, Hermanns H, Katoen JP. 2001, Berlin, Germany: Springer, 137.Google Scholar
 Kurtz TG: The Relationship between Stochastic and Deterministic Models for Chemical Reactions. J Chem Phys. 1972, 57 (7): 29762978. 10.1063/1.1678692.View ArticleGoogle Scholar
 Heiner M, Koch I, Will J: Model validation of biological pathways using Petri nets demonstrated for apoptosis. BioSystems. 2004, 75: 1028.View ArticleGoogle Scholar
 Heiner M, Mahulea C, Silva M: On the Importance of the Deadlock Trap Property for Monotonic Liveness. Proceedings of the International Workshop on Biological Processes and Petri Nets (BioPPN). 2010Google Scholar
 Murata T: Petri nets: properties, analysis, and applications. Proceedings of the IEEE. 1989, 77 (4): 541580. 10.1109/5.24143.View ArticleGoogle Scholar
 Burch JR, Clarke EM, McMillan KL, Dill DL, Hwang J: Symbolic model checking: 10^{2}0 states and beyond. Proceedings of the 5th Annual IEEE Symposium on Logic in Computer Science (LICS'90)). 1990Google Scholar
 Feller W: An Introduction to Probability Theory and its Applications, Vol. 1. 1968, John WileyGoogle Scholar
 Gillespie D: A rigorous derivation of the master chemical equation. Physica. 1992, 188: 404425. 10.1016/03784371(92)90283V.View ArticleGoogle Scholar
 Gillespie D: Exact Stochastic Simulation of Coupled Chemical Reactions. The Journal of Physical Chemistry. 1977, 81: 23402361. 10.1021/j100540a008.View ArticleGoogle Scholar
 Gillespie D: Approximate accelerated stochastic simulation of chemically reacting systems. J Comp Phys. 2001, 188: 17161733.Google Scholar
 Voit EO: Computational Analysis of Biochemical Systems. 2000, Cambridge University PressGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.