 Research
 Open Access
Analysis of intercellular signal transduction in the tumor microenvironment
 Haijun Gong^{1}Email author
https://doi.org/10.1186/175205097S3S5
© Gong; licensee BioMed Central Ltd. 2013
 Published: 16 October 2013
Abstract
Background
Recent cancer studies revealed, the interaction between pancreatic cancer cells and pancreatic stellate cells is of importance in the cancer progression. The activation of stellate cells is mediated by some growth factors and cytokines secreted by the cancer cells. In turn, the activated stellate cells will synthesize and secrete multiple growth factors to continuously stimulate the growth of surrounding cancer cells through paracrine pathways. The mechanism behind the evolution of stellate cells from quiescent state to a cancerassociated phenotype is still not well understood.
Results
To systematically investigate the interaction between cancer cells and stellate cells, we constructed a multicellular discrete value model, which is composed of several intracellular and intercellular signaling pathways that are frequently mutated in the pancreatic cancer, to study the cell cycle progression and angiogenesis. We, then, introduced and applied a formal verification technique, Symbolic Model Checking, to automatically analyze the cells' proliferation, angiogenesis and apoptosis in the proposed signal transduction model of tumor microenvironment.
Conclusions
Our studies predicted some important temporal logic properties and dynamic behaviors in the pancreatic cancer cells and stellate cells. The verification technique identified several signaling components, including the RAS, RAGE, AKT, IKK, DVL, RB and PTEN, whose mutation or loss of function can promote cell growth and inhibit apoptosis, some of which have been confirmed by existing experiments. Our formal studies demonstrated that, the bidirectional interaction between cancer cells and stellate cells could significantly increase cell proliferation, inhibit apoptosis, induce tumor angiogenesis, and promote cancer metastasis.
Keywords
 Tumor microenvironment
 signaling pathway
 discrete value model
 model checking
 formal verification
Background
Pancreatic ductal adenocarcinoma (PDAC) is a form of cancer in the pancreatic duct, which is the fourth leading cause of cancer death in the United States, and it has an extremely poor prognosis. The pathological study of PDAC has revealed a number of genetic mutations [1], including the KRAS, CDKN2A, and TP53 genes. A recent global genomic analysis [2] has identified around ten cellular signaling pathways that are frequently altered in pancreatic cancers, including the pathways of Hedgehog, Wnt, Notch, KRAS, apoptosis, TGFβ, cJUN, and G1/S phase transition. In addition, a number of growth factors and cytokines, for example, the Insulinlike growth factor (IGF), Insulin, Hedgehog (Hh), transforming growth factor (TGFβ), and the Advanced Glycation End products (AGEs) are overexpressed in the microenvironment of pancreatic cancer cells, leading to uncontrolled cancer cell proliferation, unorganized angiogenesis and evasion of apoptosis.
Recent experimental studies in pancreatic cancer [3–5] revealed, the interaction between pancreatic cancer cells (PCCs) and pancreatic stellate cells (PSCs, stromal cells of the pancreas) can stimulate cancer progression and tumor angiogenesis (formation of new blood vessels). Pancreatic cancer cells can recruit and activate PSCs to produce and maintain a growthpermissive environment for cancer progression and drug resistance. The activation of PSCs is mediated by several growth factors and cytokines, and many of which are secreted by the pancreatic cancer cells. In turn, the activated PSCs will synthesize and secrete multiple cytokines and growth factors, including Hedgehog and Wnt, through the paracrine and autocrine feedback loops to continuously stimulate cancer cells' growth. These bidirectional interactions [4] will promote cancer progression and unorganized angiogenesis. Besides, PSCs can also secrete a large amount of extracellular matrix (ECM) proteins, which are important components of the fibrous tissue along with stromal cells. Thus, the tumor microenvironment of pancreatic cancer cells and the bidirectional interaction with stellate cells can significantly increase cell proliferation, inhibit apoptosis, induce tumor angiogenesis, and promote cancer metastasis. The mechanism behind the evolution of PSCs from quiescent state to a cancerassociated myofibroblastlike phenotype is still not very clear. Several findings [4, 5] have indicated that the proangiogenic factor VEGF is of considerable importance in the stellate cell's activation and angiogenesis. To systematically understand the tumor microenvironment and the bidirectional interaction between cancer cells and stellate cells, it is imperative to investigate the intracellular and intercellular signaling pathways that regulate the cell cycle progression and angiogenesis.
Our previous work [6–9] developed Statistical Model Checking and Symbolic Model Checking techniques to study the intracellular signaling pathways in a single cell. Since the pathways implicated in the tumor microenvironment are highly interconnected, to the best of the author's knowledge, no computational multicellular model has been developed to study the interaction between pancreatic cancer cells and stellate cells due to the complexity of networks. In this work, we construct a novel in silico discrete value model of multicellular signaling pathways, which are frequently mutated [2] in pancreatic cancers, to study the interaction between PSCs and PCCs. Our 3cell model is composed of two types of cells: two pancreatic cancer cells (PCCs) and one stellate cell (PSC), which are regulated by the Hedgehog, Wnt, AGE, RbE2F, P53, RAS, PI3K, VEGF and IGF signaling pathways. Since the mechanism behind the interaction between PCCs and PSCs is not well understood, our model and analysis will provide some insights into the study of tumor microenvironment and the evolution of stellate cell from a quiescent state to an active state.
In order to formally and automatically analyze the complex network, we introduce a powerful verification technique, called Model Checking [10], which determines whether or not a model (statetransition system) satisfies a desired property expressed in a temporal logic formula. Let M be a statetransition system or a model, ${S}_{0}$ be a set of starting states, and $\psi $ be a temporal logic formula. The Model Checking problem is to verify that, for all states $s\in {S}_{0}$, the model M satisfies $\psi $  denoted by $M,s\models \psi $. Model Checker performs an exhaustive search of the state space of the model to verify or falsify the proposed temporal logic formula. Model Checking has been successfully applied to verify hardware systems and digital circuits design. In this work, finally, we apply the Symbolic Model Checking technique to analyze the complex intercellular network of pancreatic cancer cells and stellate cells. A number of important temporal logic and dynamic properties, which specify certain behaviors of regulatory components abstracted from the in vitro or in vivo experiments in the literature, are proposed to investigate the multicellular signaling pathways in the tumor microenvironment.
Methods
Multicellular model of signaling pathways
Intracellular signaling pathways
The paracrine Hedgehog (Hh) signaling is critical for the development of epithelial cells [1, 2]. In particular, Hh ligands secreted by the epithelial tumor cells can activate Hh signal transduction in the surrounding stromal cells to stimulates the cell proliferation and contributes to tumorigenesis.
Hedgehog pathway: Hh $\u22a3$ PTCH $\u22a3$ SMO $\to $ GLI $\to $ {Hh, CyclinD, ...}. The Hedgehog (Hh) ligand and its receptor Smoothened (SMO) are continuously activated or overexpressed in laterstage pancreatic carcinomas [11], while tumor suppressor protein patched (PTCH) is frequently mutated or lossoffunction, leading to a constitutive activation of Hh pathway. In the quiescent cell without Hh, SMO's activity is inhibited by forming a complex with PTCH. Once Hh binds to PTCH, SMO will be released to activate the GLI (gliomaassociated oncogene homologue) to be an active form of transcription factor. The Hh signaling pathway alone is sufficient to drive pancreatic neoplasia [12], and it is known that the activation of the HhGLI pathway is associated with tumor proliferation and pancreatic cancerrelated fibroblasts [13].
Wnt signaling pathway regulates the processes of angiogenesis and inflammation, and several proteins are genetically altered in most of pancreatic cancers according to the global genomic analysis [2].
Wnt pathway: Wnt $\to $ FZD $\to $ DVL $\u22a3$ GSK3β $\u22a3$ βCatenin $\to $ TCF $\to $ {HIF1, CyclinD, ...}. The canonical WNT pathway is activated by the interaction of Wnt and Frizzled (FZD), leading to the disassembly of AxinAPCGSK3β complex. Later, the βCatenin is translocated to the nucleus to activate the TCFLEF transcription factors [14], promoting the transcription of Cyclin D and HIF1. However, when the Wnt ligand is absent, βCatenin is localized in the cytoplasm whose activity will be inhibited by forming a complex with the Axin, APC, and GSK3β [15]. The lossoffunction or continuous activation of some regulatory components in Wnt pathway [16] is responsible for the abnormal vascular development and unorganized angiogenesis.
Recent pancreatic cancer study [17] revealed, the overexpression of the Advanced Glycation End products (AGEs), for example, HMGB1 and its receptor RAGE, is associated with the pancreatic cancer cell's survival. Our previous stochastic and deterministic simulations predicted a dosedependent p53 and Cyclin E response curve to increasing HMGB1 stimulus in a single cancer cell [6]. AGE pathway regulates the processes of inflammation and angiogenesis.
AGERASNFκB pathway: (1) AGE $\to $ RAGE $\to $ NFκBpathway; (2) AGE $\to $ RAGE $\to $ RASERKpathway. The Advanced Glycation End products (AGE), e.g., HMGB1, released by stressed or dying cells, can activate two key signaling pathways [6, 7], the RAS pathway: RAS $\to $ RAF $\to $ MEK $\to $ ERK $\to $ CyclinD, which regulates the cell cycle progression through G1 phase; and the NFκB pathway: IKK $\u22a3$ IκB $\u22a3$ NFκB $\to $ {IGF, HIF1, Hh, Wnt, AGE, ...}. In the resting cell, NFκB is located in the cytoplasm, bound to and inhibited by the tumor suppressor IκB. Once activated by AGE, the IκB kinase (IKK) will phosphorylate and deactivate IκB, leading to the translocation of NFκB into the nucleus to promote the transcription of several genes, including Cyclin D, its inhibitors A20 and IκB (frequently mutated or loss of function) [18, 19], and AGEs [20]. The overexpression of NFκB can also stimulate the secretion of VEGF through activating the HIF1 pathway.
The cell cycle progression is strictly regulated by tens of signaling pathways, and one of the hallmarks is the G1S phase transition regulated by the RBE2FCyclin E pathway. The global genomic analysis [2] identified several frequently altered regulatory components in this pathway, for example, INK4a and ARF (encoded by CDKN2A) mutations occur in 90% of pancreatic cancers [1].
G1/S phase transition pathway: {ERK, TCF, GLI, ...} $\to $ CyclinD $\u22a3$ RB $\u22a3$ E2F $\to $ CyclinE. Some upstream components of the signaling pathway, for example, ERK, TCF and GLI, can activate Cyclin DCDK4/6 complex which regulates cell cycle progression. In the normal cell, the unphosphorylated RB (a tumor suppressor protein) inhibits E2F's transcription activity by forming RBE2F complexes. E2F will be activated once its inhibitor RB is phosphorylated and inhibited by Cyclin D. E2F can promote the transcription of Cyclin E [21] and CDK2 complex, which regulates the cell cycle progression from G1 to S phase.
Drug resistance presents a challenge to the treatment of pancreatic cancer. Tumor microenvironment and angiogenesis are two important factors contributing to the drug resistance and cancer development. The environment surrounding a solid tumor is often hypoxic (low oxygen), so that angiogenesis is necessary to provide oxygen and nutrition to support the tumor's growth. The vascular endothelial growth factor (VEGF), a proangiogenetic factor whose secretion is mediated by the HIF1 pathway, can induce angiogenesis. In hypoxic conditions, HIF1 will be activated and stabilized to regulate the transcription of VEGF [14]. Moreover, it has been reported that the Wnt and NFκB pathways could also upregulate the expression of HIF1 and VEGF in the cancer cell.
Intercellualr paracrine signaling pathways
New blood vessels formation (angiogenesis) is one of the key processes in the pancreatic cancer metastasis, and this process is regulated by several proangiogenic factors, for example, VEGF, PDGF and HGF. Xu et al's study [5] revealed that the interaction between pancreatic cancer cells (PCCs) and pancreatic stellate cells (PSCs) can stimulate cancer progression and angiogenesis. The activation of PSCs could be mediated by IGF (insulinlike growth factor) and VEGF, which are produced and secreted by pancreatic cancer cells. In turn, the activated PSCs will synthesize and secrete multiple cytokines and growth factors, including Hh, Wnt, AGE, VEGF, etc., to stimulate the growth of neighboring cancer cells and promote angiogenesis through paracrine signaling pathways. In order to systematically investigate the tumor microenvironment, a simple model, which is composed of two PCCs and one PSC, was constructed to investigate the intracellular and intercellular signaling pathways that regulate the cell cycle progression and angiogenesis. In our model, the cells share similar intracellular signaling pathways that were discussed in the last section, and the bidirectional interactions is mediated by VEGF, IGF, WNT, AGE and Hedgehog pathways.
VEGFPI3KNFκB pathway: (1) VEGF $\to $ VEGFR $\to $ PI3K $\to $ PIP3 $\to $ AKT $\to $ MDM2 $\to $ P53 $\to $ Apoptosis; (2) VEGF $\to $ NFκBpathway $\to $ {Hh, Wnt, AGE, HIF1, IGF, ...}. VEGF secreted by the cancer cells can bind to its receptor VEGFR on the surrounding stellate cells to activate the PI3K pathway, which will initiate a series of reactions including the phosphorylation of PIP2, AKT, and MDM2, and repression of P53's transcription activity in the nucleus [22]. The tumorsuppressor protein P53, also called the "guardian of the genome", is mutated in more than 50% of pancreatic adenocarcinomas [1]. It is known that P53 can activate the transcription of oncoprotein MDM2 and tumorsuppressor protein PTEN, which is an inhibitor of the AKT pathway and can induce cell cycle arrest. VEGF can also activate the NFκB pathway to promote the transcription and secretion of Hh, Wnt, AGE, HIF1, IGF and VEGF, stimulating the growth of surrounding cancer cells through paracrine feedback loops.
Insulin or Insulinlike growth factor (IGF) pathway can stimulate the growth of pancreatic cancer cells and stellate cells, and inhibit apoptosis through binding and activating its receptor (IGF1R).
IGFRAS pathway: Insulin/IGF $\to $ IGFR $\to $ RAS $\to $ RAF $\to $ MEK $\to $ ERK $\to $ CyclinD. The overexpressed growth factors, e.g., the Insulinlike growth factor (IGF) and/or Insulin, can activate the RAS protein (KRAS has a high mutation frequency in pancreatic cancers [1]), leading to the phosphorylation and activation of its downstream proteins RAF, MEK, and ERK [23]. Active ERKs enter the nucleus to phosphorylate the transcription factors myc and promote the expression of the cell cycle regulatory protein Cyclin D, enabling the cell cycle progression through the G1 phase.
Discrete value model
The interaction between the pancreatic cancer cells and stellate cells is regulated by tens of proteins and crosstalk of different signaling pathways. The traditional computational techniques, including the ordinary differential equation and stochastic simulation methods [6, 7], need calculate the reaction rate of each biochemical reaction in the signal transduction. But many parameters are unknown or difficult to be estimated from existing experiments. Our aim is to qualitatively investigate the bidirectional interaction between PCC and PSC in the tumor microenvironment and compare with the experiment. In this work, we develop a discrete value model to describe the expression levels of different signaling components and dynamics of the signaling pathways without introducing any unknown parameters in the biochemical reactions.
In a discrete value model, each node represents a protein or a lipid involved in the signaling pathway. The expression level (state value) of each regulatory component (node) in the pathway can take discrete values at any specific time, for example {0, 1, 2, ..., n}, namely, 0 = "not expressed", ..., n = "overexpressed". Boolean Network is a special case of discrete value model, which can only take a Boolean value of either ON (1) or OFF (0). The evolution (state update) of each node from time t to t+ 1 is described by a discrete state transfer function, which is dependent on the state of the neighboring nodes. In this work, we assume each node can take a value of {0, 1, 2} (it can be extended to n discrete values). Similar to our previous work [8, 9, 24], the neighboring nodes are classified as activators or inhibitors: an activator node can promote or activate the expression of its downstream nodes, while the inhibitor node will inhibit or repress the expression of its downstream nodes. Since this work is the first attempt to investigate the tumor microenvironment using a computational method, for simplicity, in our discrete value model, we assume all the nodes' states are updated synchronously, i.e., the state of each node evolves according to its transfer function at any time step. This assumption has worked well in others' and our previous works on Boolean modeling [8, 9, 25, 26]. Since the biochemical processes in the cells may evolve at different rates, and the synchronous model cannot capture all the information in the cells, we will develop an asynchronous model to investigate the signal transduction in our future work.
where, ${a}_{i}$ and ${b}_{j}$ are the values of the activators ${A}_{i}$ and inhibitors ${B}_{j}$ of the node ${X}_{n}$, respectively. The values "0", "1", and "2" denote the state of "inhibited", "active", and "overexpressed" respectively. For example, PIP3 is activated by PI3K but inhibited by PTEN. At some time step t, if PI3K(t) = 1, PTEN(t) = 0, then, at the next time step, PIP3(t + 1) = 1; if PI3K(t) = 2, PTEN(t) = 0, then, PIP3(t + 1) = 2; if PI3K(t) = 1, PTEN(t) = 1, then, PIP3(t + 1) = PIP3(t).
For example, the proteins RB and E2F, under normal conditions, RB represses E2F's transcription activity by forming RBE2F complexes. When RB is phosphorylated by Cyclin D, E2F will be activated. So, RB(t + 1) = 0, if CyclinD(t) ≥ 1; else, RB(t + 1) = {1, 2}. This assumption is similar to previous work [8, 9, 25], and consistent with some clinical observation and experimental studies. Many tumor suppressor proteins, including P53, PTEN, INK4a, and ARF, are either mutated or lost in the early or late stages of PDAC. So, they cannot inhibit their downstream oncoproteins; while the oncoproteins, e.g., KRAS and NFκB, are continuously activated or overexpressed, leading to uncontrolled cell growth.
However, the output signals, including, Proliferate(a, b, c), Apoptosis(a, b, c) and Angiogenesis(a, b, c), are Boolean variables, which are activated by Cyclin E, P53 and VEGF respectively. In this work, if the corresponding activator's value is greater or equal to 1, the output signal will take a Boolean value of 1 (True); else, it will take a value of 0 (False).
The discrete value model in Figure 1 describes the interactions of different signaling components in the tumor microenvironment, which is composed of m = 96 nodes, including 3 control nodes (input signals: Hh, AGE, Wnt), and 9 output nodes (Angiogenesis, Proliferation, Apoptosis in 3 different cells). The structure depicted in Figure 1 represents a "circuit layout" of the cancerstellate cells model instead of a "state transition" system. Each node in the model is a variable taking n possible discrete values (n = 2 for 9 output boolean nodes and n = 3 for the other 87 nodes in this model), so the number of possible configurations is ${n}^{m}$ (${3}^{87}{2}^{9}$, around $1.6\times 1{0}^{44}$ possible states). When n and m are large, the network will have an astronomical amount of possible states. So, it is not realistic to use traditional computational techniques, for example, BooleaNet method [27] and stochastic simulation algorithm [6, 7], to analyze such a large network in a fast and effective way. Given a large crosstalk model of signaling pathways, one of our interests is to discover and identify some key cellular components and signal transduction sequences that will drive the system to a prespecified state (e.g., apoptosis, proliferation or angiogenesis) [8, 9, 28] at or before a prespecified time point.
We propose to apply this multicellular computational model to investigate the cellcell interactions of cancer cells with their surrounding microenvironment, in particular, with stellate cells; analyze the paracrine signaling pathways regulating the angiogenesis; identify important proteins that will drive different cells to the "apoptosis", "proliferation" and "angiogenesis" states; simulate the temporal and dynamic behaviors of the cancer cells and stellate cells in various conditions (i.e., drug interference, knockout, or overexpression). To answer these questions, we will introduce the Model Checking and temporal logic properties in the next section.
Model Checking
Model Checking is a powerful and automatic formal verification technique for finite state transition systems modeled by a Kripke structure [10], which is written as a tuple $M=\left(S,{S}_{0},R,L\right)$, where, ${S}_{0}\subseteq S$ is a set of initial states, $R\subseteq S\times S$ is a transition relation between states, and $L:S\to {2}^{AP}$ is a function that labels each state s with the set of atomic propositions (AP = {p, q, r, ...}) true in s. Given a Kripke structure M and a temporal logic formula $\psi $ expressing some desired property, the Model Checking problem [10] is to find the set of all states in S that satisfy $\psi $, i.e. to compute the set ${S}_{\psi}=\left\{s\in SM,s\models \psi \right\}$. The model M satisfies $\psi $ if ${S}_{0}\subseteq {S}_{\psi}$, otherwise, the model checker will output a counterexample (i.e., a sequence of transitions which starts from a state in ${S}_{0}$) that falsifies the formula $\psi $.
In the model checking, Computation Tree Logic (CTL) is developed to describe the properties of computation trees. The root of the computation tree corresponds to an initial state and the other nodes on the tree correspond to all possible sequences of state transitions (paths) from the root [10]. A CTL formula is constructed from atomic propositions (AP), Boolean logic connectives (i.e.,  (or, $\vee $), & (and, $\wedge $), ! (not, ¬), → (implies)), temporal operators and path quantifiers. In the CTL formula, four important temporal operators are used to describe properties on a path [6–9]: X p  p will be true in the neX t state of the path; F p  p will be true at some state in the F uture (eventually) on the path; G p  p is G lobally true (always, at every state on the path); p U q  p holds U ntil q holds. In a CTL formula, the operators X, F, G, and U must be immediately preceded by a path quantifier A  for A ll paths, or E  there E xists a path. Previous work [10] has shown that any CTL formula can be expressed in terms of ¬, $\vee $, EX, EU and EG. In this work, we proposed CTL formulas to describe the behaviors or properties of some regulatory components in the signaling pathway. For example, the formula AG($MDM{2}_{active}$ → AX$P5{3}_{inhibited}$) means, whenever an MDM2 activation event occurs, it will always inhibit P53's transcription activity in the next time step.
A path π in a Kripke structure $M=\left(S,{S}_{0},R,L\right)$ is an infinite sequence of states, that is, $\pi ={s}_{0},{s}_{1}$, ..., where ${S}_{0}$ is an initial state, ${s}_{i}\in S$ and for every i ≥ 0, $\left({s}_{i},{s}_{i+1}\right)\in R$. We use ${\pi}^{i}$ to denote the suffix of $\pi $ starting at state ${s}_{i}$. The semantics of CTL is defined as:

$M,s\models p$ iff $p\in L\left(s\right)$;

$M,s\models \neg \psi $ iff $M,s\models \psi $ does not hold;

$M,s\models {\psi}_{1}\vee {\psi}_{2}$ iff $M,s\models {\psi}_{1}$ or $M,s\models {\psi}_{2}$;

$M,s\models X\psi $ iff $M,{\pi}^{1}\models \psi $;

$M,s\models F\psi $ iff there exists a $k\ge 0$ such that $M,{\pi}^{k}\models \psi $;

$M,s\models G\psi $ iff for all $k\ge 0$, $M,{\pi}^{k}\models \psi $;

$M,s\models {\psi}_{1}\cup {\psi}_{2}$ iff there exists a $k\ge 0$ such that $M,{\pi}^{k}\models {\psi}_{2}$ and for all $0\le j<k$, $M,{\pi}^{j}\models {\psi}_{1}$;

$M,s\models E\mathrm{\Phi}$ iff there exists a path $\pi $ from s such that $M,\pi \models \mathrm{\Phi}$;

$M,s\models A\mathrm{\Phi}$ iff for every path $\pi $ from s, $M,\pi \models \mathrm{\Phi}$;
where $M,\pi \models \mathrm{\Phi}$ means, the path π in M satisfies the path formula $\mathrm{\Phi}$. More details of CTL semantics can be found in [10].
Symbolic Model Checking
Model Checking algorithm can automatically and exhaustively search the state transition system to determine, whether or not, a given model M satisfies a desired temporal logic formula $\mathrm{\Phi}$. The original Model Checking algorithm [10] represents the state transitions explicitly. It verifies or falsifies a CTL formula $\mathrm{\Phi}$ by recursively labeling the state graph with the subformulas of $\mathrm{\Phi}$, and then the graph is parsed to compute its truth value in a state for each subformula according to the CTL operators and the truth values of its subformulas [10, 24]. This algorithm could lead to a state explosion problem.
During model checking, a model or state transition system (e.g., a signaling pathway) can be described using the SMV language, and a desired cellular behavior or phenomenon can be translated into a CTL formula. Then, SMV model checker will automatically verify or falsify the CTL formula of this model. The output of the verification could be either "true" (property is satisfied) or a counterexample trace showing why the property is false (not satisfied). The complexity of the Symbolic Model Checking algorithm is O($\Psi $ (S + R)), where $\Psi $ is the size of the CTL formula, S and R are the number of states and transitions respectively [10].
Results and discussion
The multicellular model of tumor microenvironment illustrated in Figure 1 is composed of two pancreatic cancer cells (PCCs, cell A & C) and one stellate cell (PSC, cell B). The suffix a, b, c in each node represent the cell (A, B, C) that a molecule belongs to. For example, "PI3Ka" represents a "PI3K molecule" in the cancer "cell A". In this model, we use "Proliferate", "Apoptosis" and "Angiogenesis" to represent the fates of the cells, whose initial values are set to be FALSE (0). While, the other nodes, initially, can take a value of either 0 (inhibited) or 1 (active) (a stochastic effect is introduced). We apply SMV model checker to exhaustively and automatically analyze some temporal and dynamic behaviors in the signaling pathways. For simplicity, in this work, we sometimes put (a, b, c) behind the temporal logic formulas to represent some regulatory components in the cell A, B and C respectively.
In the cancer studies, we expect to predict the cancer cell's fate in various conditions (i.e., gene knockout, drug interference, loss of function, overexpression), identify key signaling components which play an important role in the tumorigenesis, and explore temporal and dynamic properties in the tumor microenvironment, in a fast and effective way. We first investigate the fates of pancreatic cancer cells (A and C) and stellate cell (B) in a predefined initial condition, that is, all the growth factors are overexpressed in the beginning.
Cell fate
Property 1: AF (Proliferate(a, b, c)); and AF (VEGF(a, b, c) ≥ 1).
Property 1 means, when the proteins Hh, Wnt, AGE (input signals) surrounding the cancer cell A are all overexpressed (i.e., initially HH(a) = 2, WNT(a) = 2, and AGE(a) = 2), the pancreatic cancer cells and stellate cell will finally reach the "Proliferate" state, and VEGF molecules in all three cells (A, B, C) will be in a state of "active" or "overexpressed", for all paths. These two properties are verified to be true. That is, overexpression of some growth factors will stimulate the synthesis and secretion of VEGF from the cancer cell (A), which might activate the stellate cell (B) by the paracrine signaling pathways, promoting the growth of stellate cell B and cancer cell C finally.
Property 2: a) AF (! Apoptosis(a, b, c) & P53(a, b, c) < 1); b) AF (Apoptosis(a, b, c)).
Property 2 a) tests, for all paths, whether or not, the cancer cells (A, C) and stellate cell (B) will finally reach "Apoptosis" state when Hh, Wnt, and AGE are overexpressed. SMV model checker verified that, all cells cannot reach "Apoptosis" state and the tumor suppressor P53's expression is suppressed (taking a value of 0); while the property b) is falsified (a counterexample trace is provided by SMV). Property 1 and 2 can be summarized as one property which is verified to be true using the SMV model checker:
Property 3: AF (P53 < 1 & ! Apoptosis & Angiogenesis & Proliferate)(a, b, c).
This property means, the cancer cell or stellate cell will reach a state in which apoptosis is inhibited while cell proliferation and angiogenesis are activated. This property is consistent with [5]'s discovery, and explains why the overexpressed growth factors and bidirectional interaction between pancreatic cancer cells and stellate cells can promote cancer progression and angiogenesis, and inhibit the cell's apoptosis.
Next, we will apply the SMV model checker to identify key regulatory components (the most frequently mutated genes or driver genes, etc.) and signal transduction sequences that will drive the system to a prespecified state (e.g., apoptosis, proliferation or angiogenesis) in the tumor microenvironment of pancreatic cancer. We assume the initial values for all nodes, except the output signals, can take a value of either 0 or 1.
Identification of key oncoproteins
Property 4: AG{(RAS(a) = {1, 2}) → AF(Proliferate(a))};
Property 5: AG{(RAS(a) = 2) → AF(VEGF = 2 & Proliferate & Angiogenesis & ! Apoptosis)(a, b)}.
KRAS mutation occurs in more than 90% of pancreatic cancers [1], especially in the precancerous stage. Constitutively active RAS pathway can stimulate the production of other key proteins during the tumorigenesis. Property 4 verified that, if the oncoprotein RAS in the cancer cell A is active or overexpressed (taking a value of 1 or 2), cell A will finally reach a "Proliferate" state, for all paths. If RAS is mutated or overexpressed (with a value 2) in the cancer cell A, it will stimulate the production and secretion of VEGF, which promotes both cancer cell A and stellate cell B to eventually enter the "Proliferate & Angiogenesis" state, while "Apoptosis" is inhibited. The verified "Property 5" demonstrated that, abating (or turning "off") the signaling function of RAS could provide a rational therapy for pancreatic cancer, and paracrine pathways play an important role in mediating the PSCPCC interaction in the tumor microenvironment.
Using model checker, besides RAS protein, we can propose similar properties to identify other key oncoproteins whose constitutive activation or mutation in the corresponding signaling pathways will influence the cell's fate. Several oncoproteins, including RAGE, DVL, AKT, and IKK, were verified to play an important role in the tumorigenesis. The following properties were checked.
Property 6: AG{(RAGE(a) = 2) → AF(NFκB = 2 & Proliferate & Angiogenesis & ! Apoptosis)(a, b)};
Property 7: AG{(DVL(a) = 2) → AF(NFκB = 2 & Proliferate & Angiogenesis & ! Apoptosis)(a, b)};
Property 8: AG{(AKT(a) = 2) → AF(VEGF = 2 & Proliferate & Angiogenesis & ! Apoptosis)(a, b)};
Property 9: AG{(IKK(a) = 2) → AF(HIF1 = 2 & Proliferate &Angiogenesis & ! Apoptosis)(a, b)}.
Property 6 and 7 predicted that, overexpression of RAGE or Dishevelled (DVL) will promote the expression of NFκB in both types of cells. This is consistent with Kang et al's discovery [17], expression of the receptor for advanced glycation endproducts (RAGE) can limit apoptosis and promote pancreatic cancer cell's survival. The oncoprotein AKT and IKK's expression is elevated in many cancers [32, 33]. Our previous work [6, 7], using stochastic simulation and ordinary differential equation methods, predicted a dosedependent P53, NFκB and CyclinE response curve to the increase of AKT and IKK. Property 8 and 9, using SMV model checker, revealed that, overexpression of AKT and IKK can increase the production and secretion of VEGF and HIF1 (Hypoxiainducible factor 1), promote the cancer cell and stellate cell to the "Proliferate & Angiogenesis" state and inhibit "Apoptosis". These properties suggest some possible ways to inhibit tumor growth and promote apoptosis through inhibiting AKT and IKK pathways, e.g., using the AKT kinase inhibitor (such as the drug GSK690693) and IKK inhibitor (e.g., Manumycin A).
Identification of key tumor suppressors
The cell cycle progression is regulated by both oncoproteins and tumor suppressors. Next, we apply SMV model checker to identify key tumor suppressors whose activation can promote apoptosis and inhibit proliferation. We analyzed the following properties:
Property 10: EG{(RB(a) = 2) → AF(ERK = 0 & Apoptosis & ! Proliferate & ! Angiogenesis)(a, b)};
Property 11: EG{(PTEN(a) = 2) → AF(CyclinD = 0 & Apoptosis & ! Proliferate & ! Angiogenesis) (a, b)};
Property 10': AG{(RB(a) = 2) → AF(ERK = 0 & Apoptosis & ! Proliferate & ! Angiogenesis)(a, b)};
Property 11': AG{(PTEN(a) = 2) → AF(CyclinD = 0 & Apoptosis & ! Proliferate & ! Angiogenesis) (a, b)}.
Property 10 and 11 were verified to be true, which means, in the RB or PTENtreated cells, there EXISTS a path, such that both cancer cell and stellate cell could reach the "Apoptosis" state finally, and the oncoprotein ERK and Cyclin D's expression is repressed. It explained why some singlegene targeted therapies had antitumor effects in some preclinical studies. However, property 10' and 11' were falsified by the SMV model checker, which means, targeting RB or PTEN in the cancer cell can NOT, for ALL paths, eventually promote the cells to enter a state that "Apoptosis" is ON and "Proliferate & Angiogenesis" are OFF. These properties demonstrate that, the crosstalk between different signaling pathways may be responsible for the pancreatic cancer cell survival even if some pathways are blocked by certain singlegene targeted therapies.
Necessary checkpoint
Property 12: !E{(CyclinD ≤ 1 & P53 ≥ 1 & HIF1 ≤ 1) U (! Apoptosis & Proliferate & Angiogenesis)}(a, b);
Property 12': !E{!(CyclinD >1 $\vee $ P53 < 1 $\vee $ HIF1 >1) U (!Apoptosis & Proliferate & Angiogenesis)}(a, b);
Property 13: !E{VEGF(a) ≤ 1 U (! Apoptosis & Proliferate & Angiogenesis)(b)}.
Here, we want to identify a necessary checkpoint that the pancreatic cancer and stellate cell will go through before they reach a predefined state. A possible checkpoint encoded in the Property 12 and 12' is verified to be true. It is worth to note that, property 12 and 12' are equivalent. And in this property, the operator "U" means "until". This formula means, there is no path in which the state "!Apoptosis & Proliferate & Angiogenesis" (state S2) is satisfied without satisfying "CyclinD > 1 $\vee $ P53 < 1 $\vee $ HIF1 > 1" (state S1) first. In other words, S1 is a necessary checkpoint for S2. This property demonstrated that, before reaching the cancerous state, the tumor suppressor P53 should have lost functions or been repressed, while oncoproteins Cyclin D or HIF1 are overexpressed or continuously activated in the cells. This property is consistent with existing experimental results that P53 is frequently mutated and Cyclin D is overexpressed in many pancreatic cancers [34]. Property 13 is false, which means, VEGF secreted by the cancer cell is not a checkpoint before the stellate cell reaches proliferation state.
Finally, we apply the SMV model checker to analyze some dynamic behaviors in the multicellular network. Oscillation is an interesting phenomenon in the signaling pathway, which has been studied in the single cell models [6–9] due to the existence of negative feedback loops.
Dynamic behaviors
Property 14: AG{(P53 > 1 → AF(MDM2 > 1)) &(MDM2 > 1 → AF(P53 < 1))}(a, b, c);
Property 15: AG{HH(a) > 1 → AG((P53 > 1 → AF(MDM2 > 1)) &(MDM2 > 1 → AF(P53 < 1)))(a, b, c)};
Property 16: AG{AGE(a) > 1 → AG((P53 > 1→AF(MDM2 > 1)) &(MDM2 > 1 → AF(P53 < 1)))(a, b, c)};
Property 17: AG{WNT(a) > 1 → AG((P53 > 1 → AF(MDM2 > 1)) &(MDM2 > 1 → AF(P53 < 1)))(a, b, c)};
Property 18: AG{(HH > 1 $\vee $ WNT > 1 $\vee $ AGE > 1)(a) → AG((P53 > 1 → AF(MDM2 ≥ 1)) & (MDM2 > 1 → AF(P53 < 1)))(a, b, c)}.
Recent experimental study [35] in a single cell observed a dynamic phenomenon of P53 and MDM2, whose expression levels in the nucleus continuously oscillated for more than 72 hours following γ irradiation. This phenomenon was studied in our previous statistical model checking based on stochastic simulations [6, 7] and Boolean network models in a single cell in response to HMGB1 stimulus [8]. Property 14 demonstrates that, this phenomenon also exists in the discrete value model of cancer cells and stellate cells due to a selfcontained negative feedback loop. Moreover, our multicellular model predicts that (Properties 1518), the external stimulus, for example, overexpression of Wnt, Hedgehog and AGE molecules around the cancer cell, can also induce the oscillation of P53 and MDM2's expression levels in the nucleus in the surrounding stellate and cancer cells. Properties 1518 were verified by the SMV model checker. Compared with [6, 7], the oscillation phenomenon is parameterindependent in our discrete value model using the Symbolic Model Checking method.
Conclusions
In this work, we developed a discrete value model of multicellular signaling pathways to study the interactions between pancreatic cancer cells and pancreatic stellate cells. The model incorporates several signaling pathways that are frequently mutated in the pancreatic cancer. The powerful Symbolic Model Checking technique is introduced and applied to analyze and validate this model formally. Several interesting temporal logic properties, which encode the cell fate, proteinprotein interaction and dynamic behaviors of some regulatory components, are proposed and verified. Compared with our previous statistical model checking work based on stochastic simulations [6, 7] and Boolean network method [8, 9], the beauty of this technique lies in its flexibility and universality. The signaling components in the model can take any kind of discrete values (3 possible values in this work), and it is easy to be extended to n possible values. Without introducing any unknown parameters, the proposed technique has checked up to $1{0}^{44}$ possible states of the multicellular network in tens of minutes, which is not realistic in the traditional simulation methods based on Gillespie's stochastic simulation algorithm and ordinary differential equations. Moreover, the Statistical Model Checking algorithm [6, 7] can only verify that a property is true with a probability, and it cannot output a counterexample if some property is not satisfied.
This work identified several genes or proteins, including RAS, RAGE, AKT, DVL, IKK, RB and PTEN, whose mutation or loss of function could promote the cancer cell and stellate cell's proliferation and inhibit apoptosis, leading to uncontrolled growth and unorganized angiogenesis in the future. The verified properties also explained, why certain singlegene targeted therapies, for example, the RB and PTENtreatment, might not always inhibit the growth of pancreatic cancer cells, due to the crosstalk of different signaling pathways, even if some pathway is blocked. These properties are either consistent with existing experimental studies, or could be verified or falsified by the future experiments.
We also investigated the dynamic behaviors in the PCCs and PSCs. The expression levels of P53 and MDM2 have been shown to oscillate in a single cell in the previous experimental study and our stochastic simulations [6, 7]. This work verified that, in response to external stimulus, the P53MDM2 network oscillation also exists in the discrete value model of multicellular signaling pathways. Our work revealed, the bidirectional interaction would continuously stimulate the neighboring cell's growth through activating the paracrine signaling pathways, in particular, VEGF pathway. Using Model Checking technique and discrete value model, we can only qualitatively compare with existing experimental discoveries. But formal analysis of this multicellular model still provides valuable information about the interactions between pancreatic cancer cells and stellate cells.
Since the proposed model is only composed of the signaling pathways that are frequently altered in the pancreatic cancer, we are far from capturing all the information in the tumor microenvirnoment, which is, in fact, regulated by tens of signaling pathways and hundreds of proteins. Experimental study [5] found that, the pancreatic stellate cells could secrete a large amount of extracellular matrix (ECM) proteins, which are important components of the fibrous tissue in the pancreatic cancer progression. Since this work attempts to investigate the interaction between PCCs and PSCs for the first time, we only consider the Hedgehog, WNT, AGE, VEGF and IGF proteins secreted by PCCs and PSCs, ECM was not incorporated into our model. A larger network of multicellular signal transduction in the tumor microenvironment will be explored in our future work. Moreover, in this work we assume all the reactions occur synchronously, i.e., the state of each protein (node) is updated at the same time. The synchronous model works well in this work, several interesting properties are consistent with existing experiment. However, biochemical processes may evolve at different rates, sometimes, the synchronous model cannot correctly describe the temporal and dynamic behaviors in the cell. We plan to apply Model Checking to study an asynchronous model in the future work. With the help of Model Checking, a comprehensive understanding of the signaling networks and their crosstalk will help cancer researchers to develop effective multigene targeted therapies for the pancreatic cancer patients.
Declarations
Acknowledgements
HG was supported by the new faculty startup grant from the Saint Louis University.
Declarations
Publication of this article was funded by the Saint Louis University.
This article has been published as part of BMC Systems Biology Volume 7 Supplement 3, 2013: Twelfth International Conference on Bioinformatics (InCoB2013): Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/7/S3.
Authors’ Affiliations
References
 Bardeesy N, DePinho RA: Pancreatic cancer biology and genetics. Nature Reviews Cancer. 2002, 2 (12): 897909. 10.1038/nrc949.View ArticlePubMedGoogle Scholar
 Jones S, Zhang X, Parsons D: Core signaling pathways in human pancreatic cancers revealed by Global genomic analyses. Science. 2008, 321: 18011806. 10.1126/science.1164368.PubMed CentralView ArticlePubMedGoogle Scholar
 Bachem MG, Zhou S: Pancreatic stellate cellsrole in pancreas cancer. Langenbecks Arch Surg. 2008, 393: 891900. 10.1007/s0042300802795.View ArticlePubMedGoogle Scholar
 Vonlaufen A, Joshi S: Pancreatic stellate cells: Partners in crime with pancreatic cancer cells. Cancer Res. 2008, 68: 20852093. 10.1158/00085472.CAN072477.View ArticlePubMedGoogle Scholar
 Xu Z, Vonlaufen A, Phillips P: Role of Pancreatic stellate cells in pancreatic cancer metastasis. The American Journal of Pathology. 2010, 177 (5): 25852596. 10.2353/ajpath.2010.090899.PubMed CentralView ArticlePubMedGoogle Scholar
 Gong H, Zuliani P, Komuravelli A, Faeder JR, Clarke EM: Analysis and Verification of the HMGB1 Signaling Pathway. BMC Bioinformatics. 2010, 11 (Suppl 7): S1010.1186/1471210511S7S10.PubMed CentralView ArticlePubMedGoogle Scholar
 Gong H, Zuliani P, Komuravelli A, Faeder J, Clarke E: Computational Modeling and Verification of Signaling Pathways in Cancer. Proceedings of Algebraic and Numeric Biology, LNCS. 2012, 6479:Google Scholar
 Gong H, Wang Q, Zuliani P, Lotze MT, Faeder JR, Clarke EM: Symbolic model checking of the Signaling Pathway in pancreatic cancer. Proceedings of the International Conference on Bioinformatics and Computational Biology (BICoB). 2011Google Scholar
 Gong H, Zuliani P, Clarke E: Model Checking of a DiabetesCancer Model. 3rd International Symposium on Computational Models for Life Sciences. 2011Google Scholar
 Clarke EM, Grumberg O, Peled DA: Model Checking. 1999, MIT PressGoogle Scholar
 Thayer S, Di Magliano M, Heiser P: Hedgehog is an early and late mediator of pancreatic cancer tumorigenesis. Nature. 2003, 425: 851856. 10.1038/nature02009.PubMed CentralView ArticlePubMedGoogle Scholar
 di Magliano M, Sekine S, Ermilov A: Hedgehog/Ras interactions regulate early stages of pancreatic cancer. Genes & Development. 2006, 20: 31613173. 10.1101/gad.1470806.View ArticleGoogle Scholar
 Walter K, Omura N, Hong SM, Griffith M, Vincent A, Borges M, Goggins M: Overexpression of Smoothened Activates the Sonic Hedgehog Signaling Pathway in Pancreatic CancerAssociated Fibroblasts. Clinical Cancer Research. 2010, 16 (6): 17811789. 10.1158/10780432.CCR091913.PubMed CentralView ArticlePubMedGoogle Scholar
 Vogelstein B, Kinzler K: Cancer genes and the pathways they control. Nature Medicine. 2004, 10: 789799. 10.1038/nm1087.View ArticlePubMedGoogle Scholar
 Wodarz A, Nusse R: Mechanisms of Wnt signaling in development. Annu Rev Cell Dev Biol. 1998, 14: 5988. 10.1146/annurev.cellbio.14.1.59.View ArticlePubMedGoogle Scholar
 Zeng G, Germinaro M, Micsenyi A: Aberrant Wnt/betacatenin signaling in pancreatic adenocarcinoma. Neoplasia. 2006, 8: 279289. 10.1593/neo.05607.PubMed CentralView ArticlePubMedGoogle Scholar
 Kang R, Tang D, Schapiro NE, Livesey KM: The receptor for advanced glycation end products (RAGE) sustains autophagy and limits apoptosis, promoting pancreatic tumor cell survival. Cell Death and Differentiation. 2009, 17 (4): 666676.PubMed CentralView ArticlePubMedGoogle Scholar
 van Beijnum JR, Buurman WA, Griffioen AW: Convergence and amplification of tolllike receptor (TLR) and receptor for advanced glycation end products (RAGE) signaling pathways via high mobility group B1. Angiogenesis. 2008, 11: 9199. 10.1007/s1045600890935.View ArticlePubMedGoogle Scholar
 Hoffmann A, Levchenko A, Scott M, Baltimore D: The IκBNFκ B Signaling Module: Temporal Control and Selective Gene Activation. Science. 2002, 298: 12411245. 10.1126/science.1071914.View ArticlePubMedGoogle Scholar
 Kang R, Tang D: The Receptor for Advanced Glycation Endproducts (RAGE) Protects Pancreatic Tumor Cells Against Oxidative Injury. Antioxidants and Redox Signaling. 2011, 15 (8): 21752184. 10.1089/ars.2010.3378.PubMed CentralView ArticlePubMedGoogle Scholar
 Yao G, Lee TJ, Mori S, Nevins J, You L: A bistable RbE2F switch underlies the restriction point. Nature Cell Biology. 2008, 10: 476482. 10.1038/ncb1711.View ArticlePubMedGoogle Scholar
 Haupt Y, Maya R, Kasaz A, Oren M: Mdm2 promotes the rapid degradation of p53. Nature. 1997, 387: 296299. 10.1038/387296a0.View ArticlePubMedGoogle Scholar
 Downward J: Targeting Ras Signalling Pathways in Cancer Therapy. Nature Reviews. 2002, 3: 1122.Google Scholar
 Gong H, Zuliani P, Clarke EM: Model checking of a synchronous diabetescancer logical network. Current Bioinformatics. 2013, 8: 915.Google Scholar
 Garg A, Cara AD: Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics. 2008, 24: 19171925. 10.1093/bioinformatics/btn336.PubMed CentralView ArticlePubMedGoogle Scholar
 Mendoza L, Xenarios I: A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theoretial biology and medical modeling. 2006, 3: 1310.1186/17424682313.View ArticleGoogle Scholar
 Albert I, Thakar J, Li S, Zhang R, Albert R: Boolean network simulations for life scientists. Source Code for Biology and Medicine. 2008, 3: 1610.1186/17510473316.PubMed CentralView ArticlePubMedGoogle Scholar
 Langmead CJ: Generalized Queries and Bayesian Statistical Model Checking in Dynamic Bayesian Networks: Application to Personalized Medicine. CSB. 2009, 201212.Google Scholar
 McMillan KL: PhD thesis: Symbolic model checkingan approach to the state explosion problem. 1992, Carnegie Mellon UniversityGoogle Scholar
 Bryant R: Graphbased algorithms for boolean function manipulation. IEEE Tran on Computers. 1986, 35 (8): 677691.View ArticleGoogle Scholar
 SMV Code. [http://cs.slu.edu/~gong/PSC.zip]
 Altomare D, Wang H, Skele K, Rienzo AD, KleinSzanto A, Godwin A, Testa J: AKT and mTOR phosphorylation is frequently detected in ovarian cancer and can be targeted to disrupt ovarian tumor cell growth. Oncogene. 2004, 23: 58537. 10.1038/sj.onc.1207721.View ArticlePubMedGoogle Scholar
 Eddy S, Guo S: Inducible IkBkinase/IkB kinase expression is induced by CK2 and promotes aberrant Nuclear FactorkB activation in breast cancer cells. Cancer Research. 2005, 65: 1137511383. 10.1158/00085472.CAN051602.View ArticlePubMedGoogle Scholar
 Chung DC, Brown SB, GraemeCook F: Overexpression of Cyclin D1 Occurs Frequently in Human Pancreatic Endocrine Tumors. The Journal of Clinical Endocrinology Metabolism. 2000, 85: 43734378. 10.1210/jc.85.11.4373.PubMedGoogle Scholar
 GevaZatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, Dekel E, Yarnitzky T, Liron Y, Polak P, Lahav G, Alon U: Oscillations and variability in the p53 system. Mol Sys Biol. 2006, 2: 2006.0033Google 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.