Computational Modeling of Complex Protein Activity Networks

Abstract

Because of the numerous entities interacting, the complexity of the networks that regulate cell fate makes it impossible to analyze and understand them using the human brain alone. Computational modeling is a powerful method to unravel complex systems. We recently described the development of a user-friendly computational tool, Analysis of Networks with Interactive MOdeling (ANIMO). ANIMO is a powerful tool to formalize knowledge on molecular interactions. This formalization entails giving a precise mathematical (formal) description of molecular states and of interactions between molecules. Such a model can be simulated, thereby in silico mimicking the processes that take place in the cell. In sharp contrast to classical graphical representations of molecular interaction networks, formal models allow in silico experiments and functional analysis of the dynamic behavior of the network. In addition, ANIMO was developed specifically for use by biologists who have little or no prior modeling experience. In this chapter, we guide the reader through the ANIMO workflow using osteoarthritis (OA) as a case study. WNT, IL-1β, and BMP signaling and cross talk are used as a concrete and illustrative model.

1. Introduction

1.1. Signal transduction networks

At any given point in time, cells are exposed to many different signals from their environment. Cells will have to interpret this multitude of signals they receive. Signal transduction networks relay and integrate signals from membrane-bound receptors, via protein activation, to the nucleus in order to regulate cellular processes such as gene transcription, metabolism, proliferation, differentiation, and apoptosis (programmed cell death).

Kinases play a key role in signal transduction by transferring phosphate groups to their substrates in a process called phosphorylation [1–4]. In this context, phosphorylation is basically a way to hand over a signal. In practice, kinases function by phosphorylating serine, threonine, or tyrosine residues on downstream substrates, thereby inducing conformational changes and/or charge alterations, resulting in modulation of protein activities [5].

Signal transduction pathways are connected to other signal transduction pathways in a mechanism we call cross talk. Due to this cross talk, signaling pathways are part of extensive signaling networks. Ultimately, dynamic changes in the signaling network determine cell fate. Insight into this network-regulating cell fate is important for controlling stem cell differentiation, understanding diseases such as cancer and osteoarthritis (OA), defining better diagnostics based on biomarker expression, and designing precision therapies.

1.2. Network topology and dynamics

To understand signaling networks, graphical representations are very useful (and widely used). In such graphical network representations, network topology and protein interactions are displayed in a static way. This is very useful for understanding network topology but fails to show the dynamics of the network interactions. In addition, as networks become large with many interactions between signaling molecules, it becomes harder to comprehend and predict the speed of the network interactions. Since we want to understand the dynamics of signaling networks, we need to incorporate quantitative aspects like activity levels and the timing of interactions. Understanding the interplay between the quantitative dynamics and the distributed and concurrent nature of networks with large numbers of components is a formidable task; this task can only be successfully undertaken by using methods and techniques that are adequately supported by software tools.

1.3. Computational modeling of signaling networks

The systems biology approach to understanding biological systems starts off from a scientific question and then follows an empirical cycle—or rather a positive spiral—of knowledge/theory → model → hypotheses → experiments → observations → update and/or refinement of knowledge/theory, until an answer to the original question is found (Figure 1). The model plays a pivotal role in this cycle:

1. To organize data and store knowledge

2. To structure reasoning and discussion

3. To perform in silico experiments and derive hypotheses

Figure 1.

The empirical spiral: applying the empirical cycle in successive rounds leads to a gradual buildup of knowledge.

An in silico model is always a simplified representation of biological reality and is never the aim in itself. Rather, it is a powerful means in the process of gaining an understanding of the biological system. Given its role in the empirical cycle, the process of modeling is especially effective when applied by the experts with respect to a certain biological system. Biologists usually have a good sense of cause-and-effect relationships of molecular interactions. In addition, they are the most knowledgeable on the network topology and the dynamics of the biological system they are studying. Since they also benefit most from the generation of hypotheses and from an efficient experimental design, biologists would be the primary candidates to construct models of their research topic.

As models are a formalization of knowledge or theories, an underlying formalism is needed to express this knowledge. Different formal methods have been successfully applied to construct representations of biological systems. Among these methods are Boolean logic [6, 7], ordinary differential equations (ODEs, reviewed in Ref. [8]), interacting state machines [9, 10], process calculi [11, 12], timed automata [13–15], and Petri nets [16, 17]. Most of these formal methods have been implemented into software tools to aid the process of modeling.

Mastery of most existing modeling tools requires training and experience in mathematical modeling. In this respect, a lack of tradition in quantitative reasoning and formal methods within the biological community at large is still a stumbling block for widespread application of modeling of biological systems. To overcome this, we built an intuitive method for the construction of formal in silico models of the dynamics of molecular networks, supported by a user-friendly modeling tool, (Analysis of Networks with Interactive MOdeling (ANIMO) [18]).

1.4. ANIMO

ANIMO is an activity network tool, built as a plug-in to the network visualization program Cytoscape [19] and founded on the formalism of timed automata [13–15], but does not require the user to have any previous training in formal methods [20]. This provides the advantages of formal models (in silico experiments and model checking) without renouncing to usability.

Nodes in an ANIMO network represent an activity level of any given biological entity, e.g., proteins directly involved in signal transduction (e.g., kinases, growth factors, cytokines, genes, and mRNA. An activity level is associated to each node, to represent, for example, the relative amount of phosphorylated kinase or the concentration of mRNA. The activity level of a node can be altered by interactions with other nodes. ANIMO networks can include activations (→) and inhibitions (─┤), which will increase (resp. decrease) the activity level of the target node if the source node is active. For example, A → B will increase the activity level of B if A is active. The speed at which an interaction occurs is defined by its k parameter, which can be estimated qualitatively by choosing among a predefined set of options (very slow, slow, medium, fast and very fast) or by directly inputting a numerical value. We note that using the qualitative choices already leads to useful models: choosing, for example, a slow interaction to represents the production of a protein, and a fast one for a posttranslational modification such as phosphorylation is already enough to provide a realistic behavior in a network with the proper node topology [18, 20, 21].

A finer control on the network dynamics can be obtained by choosing for each interaction an approximated scenario which allows to describe the interaction. A choice is available among three scenarios:

• Scenario 1 (default): the interaction rate is linearly dependent on the k parameter and the activity level of the upstream node. This is the simplest scenario and is advised for all interactions when first building an ANIMO model.

• Scenario 2: the interaction rate depends on the k parameter and on the activity levels of both nodes. In particular, it is linearly dependent on the activity level of the upstream node and inversely dependent on the activity of the downstream node. This scenario is used to model reactions where the availability of substrate is a limiting factor.

• Scenario 3 (AND gate): the interaction rate depends on the k parameter and on the activity level of two user-defined nodes. The user can determine whether the dependency on a node’s activity is linear or inverse. This scenario can be used to represent Boolean AND gates, such as “A AND B → C,” where it is required that both A and B are active in order for C to become active.

Additionally, the k parameter can be manually set to numerical values, expanding the default qualitative choices. Methods for parameter fitting are also present in ANIMO, which allow to automatically adapt the parameters to a given data set [22]. These features are useful when comparing a model to experimental data and allow to easily try different parameter settings before needing to extend a model with new nodes or interactions.

1.5. Experimental requirements

Biological events can often be interpreted as changes in activity. Activities could be defined as changes in concentration, phosphorylation, or localization of a protein or changes in gene expression because they are causal factors with respect to downstream effects. Therefore, the state or concentration of the molecules can be described in terms of an activity. The more active the molecule is, the stronger it will affect downstream processes.

Experimental design could be performed according to these guidelines:

1. A negative control (t = 0) is needed to determine background activity levels.

2. Including a positive control for each of the measurements indicates the potential maximum intensity in the biological system. This allows scaling of activity between 0 and 100% to construct a nondimensional model, omitting the need for precise intracellular concentrations.

3. Include measurements of molecules that either have downstream effects in the model or can be used as an output of the model.

4. Include overlapping treatment conditions to normalize experimental data between different days or assay batches.

5. Include multiple time points to provide insight into the dynamic behavior of the system. To decide how many time points should be measured and what the optimum time range is, consider the following factors. Ideally, measurements are obtained at time points starting from t = 0 until the system reaches a steady state. For most primary effects in signal transduction networks, this means measuring more time points in the first 2–30 minutes after stimulation. When peak dynamics are expected, three time points are the absolute minimum to describe each peak, one before the peak, one as close as possible to the actual peak, and one after the peak. Five time points and more allow finding and describing a peak in more detail, especially in the presence of experimental noise. If no peak dynamics are expected, at least four time points should be measured. Try to avoid having the highest measurement value as the first or last value in your time series, as it will lead to uncertainty about the actual behavior of the system.

We can discern primary (or direct) effects or higher order (or indirect) effects after treatment. Indirect effects are those in which feedback is involved. For signal transduction, the primary effect occurs in time points up to 240/480 minutes. For gene expression, primary effects typically take 4–12 hours. Higher order effects occur at different time ranges, e.g., signal transduction could occur up to 24/48 hours; for gene expression involving higher order effects, for example, in the case of cell differentiation, effects can take up to several weeks.

2. Case study: ANIMO modeling of inflammatory signals in  osteoarthritis

Many diseases are multifactorial, affected by many factors including genetic predisposition, age, trauma, sex, etc. These factors influence the network topology as well as its dynamics. This is hard to capture in static networks. To guide the reader through the ANIMO workflow, we use osteoarthritis as a case study.

Osteoarthritis (OA) is a painful, disabling disease with a high prevalence, occurring in about 15% of the population. The lifetime risk of knee OA is over 40% for man and almost 50% for women (reviewed in Ref. [23]). The lack of insight into the intricate signaling network of the cartilage has prevented the identification of highly needed disease-modifying osteoarthritic drugs (DMOADs). We aim to solve this by generating a comprehensive computational model of the signaling network in the cartilage [24]. In this chapter, we describe three important pathways in cartilage and OA development as a case study.

2.1. Osteoarthritis

Articular cartilage (AC) is a highly resilient tissue that covers the surfaces at the ends of long bones and ensures the pain-free and supple movement of our joints. The cartilage is mainly composed of one single cell type, the chondrocyte, which secretes and shapes the cartilaginous matrix that is necessary for its load-bearing properties. The biomechanical properties of the cartilage are dependent mainly on the composition, as well as the integrity of its matrix [25]. Once damaged, articular cartilage has low self-repair and regenerative capabilities eventually resulting in OA. This is due to its avascular nature, lack of innervation, and the embedding of chondrocytes in a dense matrix preventing cell migration. In addition, abnormalities in the cartilage-specific matrix cause a variety of skeletal malformation syndromes as well as adult-onset degenerative disorders such as OA.

OA is the most common form of arthritis and a leading cause of mobility-associated disability. OA is characterized by degeneration of articular cartilage, typical bone changes, and signs of inflammation, particularly in end-stage disease. The current management of OA is symptomatic, aimed at reduction of pain and at the end stage of disease total joint replacement as a successful treatment option for large joints (e.g., knee and hip) [26]. These treatments, however, do not cure the disease. There are no systemic drugs that can modify the disease process. Once the cartilage is damaged, no treatment exists that can intervene effectively, and the affected joint enters a disease continuum toward osteoarthritis.

Cartilage tissue homeostasis depends on a fine balance between catabolic (breakdown) and anabolic (buildup) processes. Homeostasis is regulated by a number of signaling pathways, including BMP and WNT signaling [27–31]. The amplitude of the signaling can be fine-tuned via antagonists in the extracellular space (reviewed in Ref. [32]). Typical catabolic pathways include the inflammatory pathways, including TNFα and IL-1β.

2.2. Osteoarthritis at the molecular level

OA is a disease caused by loss of homeostasis, resulting in altered mechanical and biochemical signals. Some of the key biochemical signals are growth factors such as WNT, IL-1β, transforming growth factor beta, bone morphogenetic proteins (BMPs), and Indian hedgehog homolog [27, 33, 34]. Mechanical stress on the extracellular matrix (ECM) plays a key role in OA development [27, 33, 34]. Any changes in this complex biological system, such as those caused by injury or aging, can disrupt cartilage homeostasis and lead to either catabolism characterized by expression of, for example, matrix metalloproteases (MMPs), and aggrecanases (ADAMs), or anabolism characterized by expression of collagen II and aggrecans [33]. As OA progresses the chondrocytes start to lose their characteristic phenotype, which in some cases results in differentiating into hypertrophic chondrocytes [35], resulting in endochondral ossification, by destroying the surrounding collagen II and replacing it with collagen X. Eventually, the hypertrophic chondrocytes will recruit osteoblasts that will proceed to form an osteophyte [35]. It is important to note that OA is not a disease that will damage the whole joint evenly. Throughout the cartilage there will be cells in different stages of differentiation, ranging from seemingly healthy cells to osteophyte forming hypertrophic chondrocytes.

The direct control of chondrogenic differentiation and hypertrophy is believed to be tightly regulated by the transcriptional activity of two main transcription factors: RUNX2, a transcription factor important for the regulation of hypertrophic differentiation, and SOX9, master transcription factor for chondrogenic development [36, 37]. The exact activity of these factors seems a key in determining the outcome of the chondrocyte phenotype.

The first steps in any computational modeling workflow are to thoroughly investigate the signal transduction pathways involved in the disease and to choose which pathways will be focused on. In this example, we will show the BMP and WNT signaling pathways for their importance in cartilage development and IL-1β as an inflammatory signal involved in osteoarthritis.

2.2.1. WNT signaling in the cartilage and osteoarthritis

The canonical Wnt pathway is crucial for cell survival and OA activation. The canonical pathway is characterized by the axin/glycogen synthase kinase-3β (GSK3-β) destruction complex which maintains low intracellular concentration of the key transcriptional regulator, β-catenin [33]. WNTs bind to the Frizzled (Fz) transmembrane receptors, resulting in the recruitment of the transmembrane protein LRP5/6. This complexation leads to the phosphorylation and dissociation of the destruction complex allowing β-catenin to accumulate and translocate to the nucleus [33]. In turn, β-catenin downregulates both collagen type 2A (COL2A) and SRY-box 9 (SOX9), leading to cell dedifferentiation and proliferation [38]. The WNT pathway can be activated by IL-1β through the phosphoinositide 3-kinase (PI3K) – protein kinase B (Akt) pathway [39, 40].

2.2.2. BMP2 signaling in osteoarthritis

BMP2 signaling is key pathway in the development of both the bone and cartilage. In endochondral bone formation, it is responsible for the clustering of the mesenchymal stem cells, the acquisition of the chondrocyte phenotype, and the final differentiation into hypertrophic chondrocytes [41, 42]. This final step is stopped in order to produce adult chondrocytes [41]. BMP2 is found in both healthy and OA adult chondrocytes [42]. BMP2 signaling occurs when BMP2 binds to its type 1 and type 2 receptors, which in turn phosphorylate mothers against decapentaplegic (SMAD) homologs 1, 5, and 8. Subsequently, SMAD 1, 5, or 8 dimers could bind to the ubiquitous SMAD-4 transcription factor [41]. BMP2 signaling can upregulate Col2a, SOX9, ColX, and MMP13 gene expression [38, 41]. Once OA is advanced, BMP2 can cause hypertrophic differentiation of chondrocytes, leading to osteophyte formation [41].

2.2.3. Interleukin 1β signaling in osteoarthritis

IL-1β is a key pro-inflammatory cytokine that drives OA progression by inducing the expression of cartilage degrading enzymes such as matrix metalloproteinases (MMPs) [43, 44].

IL-1β signals by binding to the transmembrane IL-1 receptor I (IL-1RI), leading to the activation of multiple signaling pathways. The canonical IL-1β pathway signals through NF-κB, but IL-1β can also activate the p38-MAPK and JNK-MAPK pathways. The activated receptor then assembles two signaling proteins, myeloid differentiation primary response gene 88 (MYD88) and interleukin-1 receptor-activated protein kinase 4 (IRAK4). Together, the proteins form a stable IL-1-induced first signaling module and activate IRAK1 and IRAK2, which in turn activate TRAF6, PELI 1-3, TAK1, and MEKK3 [45]. IRAK1 also activates the inhibitor of nuclear factor B kinase (IKK) complex, which is necessary for the translocation of NF-κB (nuclear factor kappa-light-chain-enhancer of activated B cells) subunits to the core. The IKK complex consists out of IKK1 and IKK2 plus the regulatory subunit NF-κB essential modifier (NEMO) and phosphorylates IkB, the inhibitor of NF-κB, leading to its degradation [45, 46]. Due to the degradation of IkB, two NF-κB subunits, p50 and p65, are released and translocate into the nucleus, where they bind to conserved DNA motifs, which exist in many IL-1β responsive genes, including the genes for IL-1B [45, 46] and IκBα [45].

For the p38 MAPK and the JNK pathways, TAK1 and MEKK3 are mainly responsible, which activate MAPK kinase kinases (MKKs) 3, 4, 6, and 7 [45]. In addition, ERK1/ERK2 is also activated as a result of the activation of the IKK complex, which influences MKK1. All three MAPK pathways influence the activation protein 1 (AP-1), affecting the DNA expression of IL-1β response genes [45].

An increasing amount of data indicates the influence of IL-1β on the degeneration of extracellular matrix in the pathology of OA [47]. In addition, we have previously shown that WNT/β-catenin inhibits IL-1β-induced MMP expression in human articular cartilage [48]. Moreover, we showed that the WNT/β-catenin-regulated transcription factor TCF4 can bind to NF-κB, thereby enhancing NF-κB activity [30].

2.3. Defining an a priori network

The aim of computational modeling is not to provide the most complete representation of all the interactions in a signaling network, but to use as many interactions as needed to provide insight into cellular mechanisms. As such, models are indeed simplified representations of the real situation: one can choose a level of abstraction depending on the available information and the research question that is asked. The level of abstraction is always a trade-off between precision and feasibility.

In the case study presented here, we do not aim to build a precise model, but aim to show how building a model enables researchers of all levels of experience to visualize, summarize, and formalize models. Generating a relatively simple model in which key interactions are shown enables researchers to test and discuss various hypotheses quickly. With the obtained insight, one can then choose to validate only those hypotheses that the researchers will expect to truly yield new information. The model is then used as a backbone for the smart design of wet lab experiments rather than the trial-and-error methods that are traditionally used in the field.

2.3.1. Defining a research question

We base our research question on our analysis of the role of the individual pathways and their possible cross talk on OA development. Although the roles of WNT and BMP in cartilage and OA development are well described, the precise interactions between WNT, BMP, and IL-1β in regulating OA and chondrocyte hypertrophy are not fully understood (reviewed in Ref. [31]).

The model can be used for rapid and reiterative queries to derive and probe hypotheses such as whether IL-1β could influence cartilage homeostasis by modulating the activity of the cartilage and bone transcription factors SOX9 and RUNX2 and their downstream targets. Similarly, the role of BMP and WNT signalings, two important pathways in cartilage development and maintenance, can be explored on their potential modulatory roles on IL-1β expression and function.

To explore this example and provide a concrete guide through the ANIMO workflow, we will first draw a priori knowledge network, then formalize this network based on literature and our own data, and then perform a few simple in silico experiments that will be validated in the wet lab.

2.3.2. Drawing a priori knowledge network

To build an a priori network, we use KEGG (www.genome.jp/kegg) and WikiPathways (wikipathways.org) to decide on the topology of the proteins in the network.

We first draw an a priori network diagram that includes some of the most important factors in the signaling pathways of interest: IL-1β (based on: WikiPathways WP2637), WNT (e.g., the WNT homepage [49]), and BMP (e.g., WP1425). We include those intracellular molecules that we can actually measure in our experiments (see Sections 1.5 and 2.6), and added dashed lines to indicate other molecules important for the pathway were omitted (see Section 2.2).

In the WNT pathway, the inhibition of the destruction pathway leads to the inhibition of destruction of β-catenin, resulting in its upregulation and nuclear translocation [32]. So the net effect of WNT binding to its receptor is the increase of β-catenin activity. To simplify the model, we omit the many steps involving the double inhibition mechanisms in the WNT pathway, resulting in the simplified path WNT → WNTR → GSK3 → β-catenin → TCF/LEF.

SOX9 regulates expression of the matrix proteins collagen 2 and aggrecan. RUNX2 regulates transcription of collagens I and X and MMP13 [50, 51]. Since the activity of SOX9 and RUNX2 is key to the switch from the cartilage to hypertrophic cartilage, we included SOX9 and RUNX2 and some of their target genes in our diagram (Figure 2, 3A).

Figure 2.

A priori knowledge network of WNT, BMP, and IL-1β pathways. IL-1β canonical (blue) and noncanonical signaling (green) showing cross talk with the transcription factors RUNX and SOX9 (both light blue), transcription factors from the WNT (red), and BMP2 (purple) pathways. Solid arrows indicate direct protein interaction, and dashed arrows indicate that intermediate protein interaction is omitted because no cross talk occurs between these proteins with other pathways. The colors indicate the canonical and noncanonical pathways corresponding to the external signals WNT, BMP2, and IL-1β.

2.4. Adding dynamics to the network

Once a priori knowledge network has been chosen, it needs to be drawn in ANIMO as a collection of nodes and interactions. For each node, a maximum number of activity levels can be chosen: unless a model is extremely large, it is safe to use 100 levels for all nodes. After providing the node with a name and an initial activity (which describes the state of the node at the start of a simulation), a description can optionally be added. Descriptions can be used as rationale for the presence of a certain node in the model, for example, a node description can contain citations to literature or references to own experimental results.

When adding an interaction between two nodes, a choice for an approximation scenario and a k parameter is necessary. The k-constants in our a priori knowledge network are not taken from literature as the strength of all interactions is assumed to be equal. This is a pragmatic decision, as many actual k-values are not described for most of the protein interactions in our network. In our initial models, we assume that there are in general two types of reactions: fast reactions for posttranslational modifications, such as phosphorylation, and slow reactions where gene transcription occurs. We therefore add reactions between nodes using these two types of reaction speed with a “Scenario 1” setting. We also add auto-inhibition to indicate inhibition as described in the literature for, e.g., receptor internalization, phosphatase activity, and, in the case of NF-κB, nuclear export as regulated by IκB. A more in-depth discussion on parameters, scenario choices, and network topology can be found in Refs. [21, 22].

Based on our experience, we expect proteins directly downstream of an activated receptor, such as p38 and JNK that are downstream of IL-1R, to be most activated by phosphorylation about 15–30 minutes after stimulation, and that the activity would be decreased to the starting situation between 60 and 240 minutes after stimulation. We therefore set the initial parameters to match these assumptions.

2.5. Testing network effects of different stimuli in silico

Testing the effect of different stimuli in a computational network is performed in small steps. During each step the parameters in the model can be updated so that the dynamics of the various nodes in the network match our knowledge.

1. Step 1. What is the normal “steady-state” (=No Input) situation of the nodes represented in the network? For example, for articular chondrocytes it is known that the transcription factor SOX9 is active and that collagen 2 and aggrecan are expressed. It is also known that the WNT and IL-1β pathways are inactive and that BMP is active at a low level [52, 53]. We therefore can adjust our starting activities to these settings. This is generated in Figure 4. We display the activities of the proteins of which we plan to measure the phosphorylation, ERK1/2, GSK3, JNK, and p38 as well as the gene expression of AXIN2 and COL2A1. Initially, COL2A is expressed, indicating SOX9 activity.

2. Step 2. In a first test of the response of cells to various stimuli, we tested the presence of WNT starting from our “steady-state” model generated as described in Step 1. Since we do not starve our cells in the experiment, BMP2 is active at a low level of 20 activity units. After WNT addition we see that GSK3 becomes activated, and the activity peaks between 30 and 60 minutes after WNT addition and then trails off around 400 minutes. We see that AXIN2 becomes present between 2 and 4 hours after WNT treatment. This is probably faster than what can be expected from a newly synthesized mRNA. At the same time, due to the inactivation of SOX9 by β-catenin [54], we observe a reduction in the COL2A1 expression around 2–3 hours after WNT addition (not shown).

3. Step 3. We then tested the presence of IL-1β starting from our “steady-state” model. Again, BMP is active at a level of 20. We now observe that within 15 minutes of IL-1β addition, the three downstream kinases ERK1/2, p38, and JNK become active. In turn, these kinases activate RUNX2 [55], thereby activating its target genes. Due to the negative feedback loop from RUNX2 to JNK and p38, the activity curve is more narrow for these proteins than it is for ERK1/ERK2 [56]. GSK3 becomes slightly active via AKT activation by IL-1β. We see no reduction of COL2A mRNA expression and a transient activation of COL1 and COLX mRNA expression.

4. Step 4. Next, we want to see the effect of dual stimulation of WNT and IL-1β when starting from a healthy situation as described under Step 1. Addition of WNT and IL-1β, in the presence of 20% BMP, decreases the activity of SOX9 and therefore causes loss of COL2A1 expression. At the same time, RUNX2 becomes active, thereby activating MMP13, COL1, and COLX, ultimately leading to a hypertrophic phenotype.

5. Step 5. The next question was what is the effect of IL-1β and WNT in the presence of high BMP activity? In our model, the presence of high BMP activity is enough to prevent the loss of COL2A expression (not shown), while at the same time leading to induction of RUNX2 activity and the corresponding COL1 and COLX expression. However, we would have expected that the high levels of WNT and IL-1β would lead to reduced SOX9 activity as seen in articles describing OA (reviewed in Ref. [31]). We therefore need to adapt our network to match the literature data.

6. Step 6. It has been described that SOX9 and β-catenin influence each other’s activity [54]. Also, RUNX2 suppresses SOX9 in bone formation [57], while SOX9 suppresses RUNX2 activity [36]. We added mutual inhibitions between SOX9 and RUNX2. This showed that when RUNX2 becomes active, it suppresses the SOX9-induced COL2 expression.

These few in silico experiments provide insight into the possible cross talk between these three pathways and their effect on SOX9 and RUNX2 activity and the possible effects of upstream signals. While asking the questions, we modified the parameters of the reactions, added inhibitory loops, and checked signaling cross talk in order to match the reaction speed in the model with our literature data or own experience.

This initial model now allows us to investigate the mechanism of IL-1β in inducing cellular hypertrophy by directly regulating SOX9 and RUNX2 activity. We hypothesize that IL-1β will increase expression of hypertrophic genes by upregulating RUNX2 activity and downregulating SOX9 activity (Figure 3).

Figure 3.

The ANIMO model built to represent the cross talk between the WNT, BMP, and IL-1β pathways. (A) The initial configuration of the model. (B) An example of simulation in ANIMO: activity levels of all nodes are shown after 120 minutes of treatment with WNT + IL-1β using color coding. The node colors are indicative of their activity at the indicated time points, with green being most active, via yellow to red, which is inactive as indicated in the figure, bottom left.

2.6. Testing hypothesis by wet lab experiments

2.6.1. Designing experiment

The outcomes of the in silico experiments are used as guideline for the experimental setup for the wet lab validation. In this case, we questioned the effect of IL-1β on WNT signaling in the presence or absence of BMP signaling in the cartilage. For this we stimulated cells with IL-1β, WNT3A, or BMP2 either alone or in combinations. Parts of the data used in this chapter are published previously [21]. The other raw data can be obtained upon request.

2.6.2. Wet lab data

After the creation of the initial model with defined nodes, the next step was to obtain the data to fit into the model. This step was carried out mainly with wet lab data complemented with literature data. It is important that the analyses of wet lab data show consistency with well-known osteoarthritic cellular responses; these analyses are done prior to inclusion of experimental data into the model.

Figure 4 shows the measured and predicted activity of the various proteins of our network. We observe attenuations in, for example, p38, JNK, and ERK phosphorylation, where WNT partially inhibited the effect of IL-β on the phosphorylation of these proteins. We have already described these data [21]. When we add BMP to the cells, in combination with IL-1β and WNT, we see that in addition to the lower activity of the proteins, there is also a delay in the time by which the maximum activity is reached. This results in a delay and a reduction of the level of gene expression of all genes tested.

Figure 4.

Comparing ANIMO’s results with wet lab data: signal transduction. The results from two versions of the model are shown: the initial version with qualitative parameters (initial model) and the one obtained with ANIMO’s automatic parameter fitting feature (fitted model). Colors are indicative of activity with green being most active and red, via yellow, inactive (see Legend).

2.6.3. Comparing wet lab data to in silico data

The wet lab data obtained from the experiment were normalized and rescaled from 0 to 100 in order to be comparable with ANIMO’s simulation results. For the complete normalization procedure, we refer to Ref. [20] and ANIMO’s manual. The resulting.csv tables, together with the model, can be found online, in the link below: http://fmt.cs.utwente.nl/tools/animo/content/models/Phosphorylation.

To evaluate the accuracy of the model, we compared its simulation results against the wet lab data concentrating at first on the signal transduction part of the network (Figure 4, first two rows). The initial match was already quite close, even if the parameters used in the model were all of qualitative nature (Table 1). The heat map graphs we show in Figure 4 were obtained in ANIMO, where wet lab data can be directly compared to simulation results.

Interaction	k^*	Qualitative parameter	Scenario	k^**
ACAN --| ACAN	0.001	Very slow	1	0.001
AKT --> GSK3	0.008	Fast	1	0.01699021
AKT --| AKT	0.004	Medium	1	0.05987754
AXIN2 --| AXIN2	0.002	Slow	1	0.00016882
Beta-catenin --> TCF/LEF	0.008	Fast	1	0.008
Beta-catenin --| beta-catenin	0.002	Slow	1	0.002
Beta-catenin --| NF-κB	0.008	Fast	1	0.03027397
Beta-catenin --| SOX9	0.016	Very fast	1	0.01824899
BMP2 --> BMPRI/BMPRII	0.016	Very fast	1	0.01769909
BMPR internalization --| BMPRI/BMPRII	0.016	Very fast	2	0.01426827
BMPRI/BMPRII --> BMPR internalization	0.001	Very slow	1	0.00067129
BMPRI/BMPRII --> SMAD1/5/8	0.008	Fast	1	0.00804472
COL1 --| COL1	0.002	Slow	1	0.002
COL2A --| COL2A	0.001	Very slow	1	0.00042541
COLX --| COLX	0.002	Slow	1	0.002
ERK1/ERK2 --> RUNX2	0.004	Medium	1	0.00417643
ERK1/ERK2 --| ERK1/ERK2	0.004	Medium	1	0.0112296
GSK3 --> beta-catenin	0.016	Very fast	1	0.01599999
GSK3 --> GSK3 phosphatase	0.008	Fast	1	0.02126481
GSK3 phosphatase --| GSK3	0.016	Very fast	2	0.03690672
ID1 --| ID1	0.002	Slow	1	0.00127961
IkB --| NF-κB	0.002	Slow	1	0.00223597
IL-1B --> IL-1RI	0.008	Fast	1	0.01240816
IL1B --| IL-1B	0.002	Slow	1	0.00060578
IL-1BR internalization --| IL-1RI	0.032	Very fast +	2	0.03583533
IL-1RI --> AKT	0.004	Medium	1	0.04684802
IL-1RI --> beta-catenin	0.008	Fast	1	0.00908979
IL-1RI --> ERK1/2	0.008	Fast	1	0.0041156
IL-1RI --> IL-1BR internalization	0.004	Medium	1	0.00782616
IL-1RI --> JNK	0.008	Fast	1	0.016
IL-1RI --> NF-κB	0.016	Very fast	1	0.04698054
IL- 1RI --> p38	0.008	Fast	1	0.04107936
JNK --> RUNX2	0.004	Medium	1	0.00414343
JNK --| JNK	0.004	Medium	1	0.00055435
MMP13 --| MMP13	0.002	Slow	1	0.002
NF-κB --> IkB	0.016	Very fast	1	0.016
NF-κB --> IL-1B	0.002	Slow	1	0.00156559
NF-κB --> MMP13	0.001	Very slow	1	0.001
p38 --> RUNX2	0.004	Medium	1	0.00399137
p38 --| p38	0.004	Medium	1	0.04260553
RUNX2 --> COL1	0.002	Slow	1	0.002
RUNX2 --> COLX	0.002	Slow	1	0.002
RUNX2 --> MMP13	0.001	Very slow	1	0.001
RUNX2 --| JNK	0.004	Medium	1	0.00606149
RUNX2 --| p38	0.004	Medium	1	0.00406615
RUNX2 --| RUNX2	0.004	Medium	1	0.01181898
RUNX2 --| SOX9	0.004	Medium	1	0.004
SMAD1/5/8 --> ID1	0.002	Slow	1	0.00270426
SMAD1/5/8 --> RUNX2	0.008	Fast	1	0.008
SMAD1/5/8 --> SMAD7	0.004	Medium	1	0.00671529
SMAD1/5/8 --> SOX9	0.004	Medium	1	0.004
SMAD7 --| SMAD1/5/8	0.008	Fast	2	0.02130379
SOX9 --> ACAN	0.001	Very slow	1	0.001
SOX9 --> COL2A	0.001	Very slow	1	0.001
SOX9 --| beta-catenin	0.001	Very slow	1	0.00028855
SOX9 --| RUNX2	0.004	Medium	1	8.7016E-05
SOX9 --| SOX9	0.002	Slow	1	0.00424899
SOX9-activating signal --> SOX9	0.016	Very fast	1	0.00255911
TCF/LEF --> AXIN2	0.001	Very slow	1	0.00021457
TCF/LEF --| TCF/LEF	0.004	Medium	1	0.004
WNT --> WNTR	0.008	Fast	1	0.02304019
WNTR --> AKT	0.016	Very fast	1	0.04626554
WNTR --> GSK3	0.008	Fast	1	0.00372738
WNTR --> WNTR internalization	0.004	Medium	1	0.00600663
WNTR internalization --| WNTR	0.016	Very fast	1	0.06497972

Table 1.

Parameters for WNT, BMP, and IL-1β signaling in the initial model (k*) and the fitted model that was optimized using experimental data (k**; see Section 2.6.4).

Another feature provided by ANIMO is parameter fitting, which allows to automatically try different parameter values for the interactions in a model, comparing the simulations with a given data set. This lets the researcher check whether the model topology can be a plausible explanation of the reference wet lab data. In case no parameter set can satisfyingly match the given data, or if only a very narrow parameter choice fits well, it is likely necessary to try a different wiring of the network model.

2.6.4. Optimizing model

Our next step was to use ANIMO’s automatic parameter fitter on the model, using the wet lab data as reference. We divided the model in two subnetworks roughly corresponding to the WNT and IL-1 pathways. The two subnetworks were independently fit, based on the wet lab data for the treatments with WNT and IL-1. Dividing the model allowed us to limit the parameter space for the automatic search, making it more rapid: it took less than 3 minutes to complete. In practice, to fit a part of the ANIMO model, we disabled the part of the network we were not focusing on. We then selected those interactions whose parameters we wanted to optimize and clicked on the “Optimize k-Values” command in ANIMO’s interface. After providing the proper file with the wet lab data, we let the tool to automatically try different k-values for the interactions, comparing the model simulations with the data and determining the fitness. Once the tool could find a better fitting set of parameters, the process would terminate, showing the resulting match for the candidate parameter set. We repeated the same procedure on both WNT and IL-1 pathways, fitting the signal transduction parts to the data. We then simulated the WNT + IL-1 treatment in the model and compared it to the data, finding it was fitting already well enough (see Figure 4, first and last row).

For BMP2 signaling we optimized SMAD activity based on phosphorylation of Western blot of SMAD1/5/8. SMAD1/5/8 was most active at 15 minutes posttreatment (data not shown).

Finally, we compared the model with the polymerase chain reaction (PCR) experimental data, repeating the fitting process only on that part of the network. The resulting model parameters can be found in Table 1 (k**), while the comparison between PCR data and ANIMO mRNA node activities is shown in Figure 5. We note that the general trends were captured in the model. Differences between activities in the model vs. the experimental data are especially visible in the longer time scales, around 24 hours. This can be expected because higher order effects, such as feedback loops, which take place on longer time scales, are not included in our model.

Figure 5.

Comparing ANIMO’s results with wet lab data: protein production. The activities of nodes representing mRNA in the ANIMO model (bottom panel) were compared against wet lab PCR data for actin, Col2a, ID1, and IL-1 (top panels).

2.6.5. Validating hypothesis

Our hypothesis that IL-1β will increase expression of hypertrophic genes by upregulating RUNX2 activity and downregulating SOX9 activity can now be tested in our optimized model. For this, we can investigate, in silico, the changes in activity of SOX9 and RUNX2 and their corresponding genes, even when no gene expression is measured in the wet lab.

Our model predicts that in the presence of IL-1β and WNT, SOX9 activity will be inhibited and RUNX2 will be activated with an initial peak activity in the first hour (Figure 6). This indicates loss of cartilage homeostasis and a slight increase in hypertrophy, which is sustained in time due to the permanent increase in RUNX2 activity (not shown). So even though the initial WNT and IL-1 signal are no longer present, an increase in RUNX2 activity results in a sustained expression of collagen 1 and collagen 10 as well as MMP13, albeit at a low level. Even high levels of BMP2 cannot prevent the loss of SOX9 activity, eventually leading to hypertrophy. These data validate our hypothesis, at least in silico.

Figure 6.

In silico validation of hypothesis: IL-1β will increase expression of hypertrophic genes, COL1, COL10, and MMP13, by upregulating RUNX2 activity and downregulating SOX9 activity.

3. Discussion and conclusion

In this chapter we provide an example of a workflow for starting computational modeling based on literature and experimental data. The aim was not to make the most comprehensive model in terms of network topology, but to understand the dynamics of the network activity in terms of signaling cross talk and corresponding downstream effects.

We chose to use the software ANIMO as a plug-in in Cytoscape as it offers a user-friendly interface in which biologists can interactively create and explore computational models of signal transduction networks. This allows to gain intellectual control over the dynamic behavior of the network that is modeled. We showed that network topologies can be constructed, modified, and enhanced with a formal description of the associated dynamic behavior. The process of modeling biological network dynamics is a prerequisite for formally comparing experimental data to a priori knowledge. ANIMO can also be used in research groups to assist in the storage and transfer of knowledge on biological networks and as a guide in discussions.

Most of the plug-ins for Cytoscape are based on static analyses, for example, they make it possible to find the hubs in a network, to cluster nodes by specific features, or to associate external data sources to the network. This allows to effectively represent large quantities of information and obtain useful insight from them, but the focus is still on the “static picture.” ANIMO concentrates instead on the network dynamics: applying an abstract representation of biochemical kinetics, it allows to represent how signaling networks evolve with time under different conditions. Graphs and node colors provide the user with useful representations of the network dynamics. Further analyses are enabled by the possibility to perform model checking on the underlying timed automata model, which can be used without the need to acquire additional training in formal methods.

ANIMO describes biological entities in the network in terms of their activity. This generalizes easily into most signal transduction processes. However, this can also be used to model any process that can be abstractly modeled as a variable activity. Examples are the inclusion of processes such as receptor internalization and phosphatase activity but also inhibition of an activated protein by proteosomal degradation or nuclear export. This flexibility helps the user in describing parts of the network for which the molecular details are unknown or of lesser importance.

In the model presented here, we show how a priori knowledge network based on three signaling pathways can be constructed and tested in silico by asking questions in small steps at a time. We then showed how experimental data of a limited number of proteins and genes, at a wide range of time points, aid to optimize topology and dynamics of the proteins and mRNAs in our network. In the next step, we can prioritize and design new experiments that can be validated in the wet lab. Seeing the role of a computational model in the empirical spiral in Figure 1, the work is never finished, but each step in the cycle aids to optimize the model and hence the molecular insight into the dynamics and topology of the cellular signal transduction network.

In ANIMO we proved our hypothesis that IL-1β will increase expression of hypertrophic genes by upregulating RUNX2 activity and downregulating SOX9 activity. For this we used a combination of literature and experimental data to optimize the model parameters. This allowed us to obtain insight into the order of events in the presence of WNT and/or IL-1β at the level of SOX9 and RUNX2 activity. In addition, it allowed insight into the complex interconnectivity of three individual pathways. Such models also yield high content data at high temporal resolution, a feat that is difficult to achieve using only wet lab approaches.

Interestingly, in one computational model, we are able to show a combination of protein activity (phosphorylation) and subsequent mRNA expression. This is a combination model of events at very different time lines. The advantage of our strategy, which included automatic parameter fitting, is the possibility to predict cell fate based on both changes in phosphorylation/protein activity and corresponding gene expression differences.

In the future, ANIMO and related tools may lead to a new paradigm for interactive representation of biological networks. Networks in digital textbooks and articles could be displayed as animations amenable to modifications by readers. Repositories of formal descriptions of signaling modules could be used to put together executable signaling networks. A more user-friendly way of interacting with dynamic network models will lead to a more thorough understanding of biological networks and will accelerate hypothesis-driven research.

Acknowledgments

We thank Patrick Kaufhold and Dr. Jacqueline Plass for the technical assistance. We thank Ricardo Urquidi Camacho for the technical assistance and valuable discussion during the time that he was an MSc graduate student in our lab in 2012 and 2013. We thank Kannan Govindaraj for critically reading the manuscript and for the valuable discussion.