- Open Access
A systems biology approach to suppress TNF-induced proinflammatory gene expressions
Cell Communication and Signaling volume 11, Article number: 84 (2013)
Tumor necrosis factor (TNF) is a widely studied cytokine (ligand) that induces proinflammatory signaling and regulates myriad cellular processes. In major illnesses, such as rheumatoid arthritis and certain cancers, the expression of TNF is elevated. Despite much progress in the field, the targeted regulation of TNF response for therapeutic benefits remains suboptimal. Here, to effectively regulate the proinflammatory response induced by TNF, a systems biology approach was adopted.
We developed a computational model to investigate the temporal activations of MAP kinase (p38), nuclear factor (NF)-κB, and the kinetics of 3 groups of genes, defined by early, intermediate and late phases, in murine embryonic fibroblast (MEF) and 3T3 cells. To identify a crucial target that suppresses, and not abolishes, proinflammatory genes, the model was tested in several in silico knock out (KO) conditions. Among the candidate molecules tested, in silico RIP1 KO effectively regulated all groups of proinflammatory genes (early, middle and late). To validate this result, we experimentally inhibited TNF signaling in MEF and 3T3 cells with RIP1 inhibitor, Necrostatin-1 (Nec-1), and investigated 10 genes (Il6, Nfkbia, Jun, Tnfaip3, Ccl7, Vcam1, Cxcl10, Mmp3, Mmp13, Enpp2) belonging to the 3 major groups of upregulated genes. As predicted by the model, all measured genes were significantly impaired.
Our results demonstrate that Nec-1 modulates TNF-induced proinflammatory response, and may potentially be used as a therapeutic target for inflammatory diseases such as rheumatoid arthritis and osteoarthritis.
The tumor necrosis factor (TNF), first termed in 1962 , was initially known for its ability to induce programmed cell death or apoptosis. As a result, throughout the years, the TNF has been intensely investigated for its anticancer property . Today, this cytokine is central to the regulation of myriad important cellular processes such as proliferation, differentiation, growth, and the immune response.
TNF binds to two types of outer membrane bound receptors on target cells, TNFR1 and TNFR2, and triggers the cell survival and proinflammatory NF-κB and MAP kinases activations . In addition, the TNFR1 induces intracellular cell death pathways via caspases after internalization through endocytosis. It is, therefore, conceivable that the dysregulation of the TNF signaling process will misbalance proinflammatory and/or apoptotic responses. Notably, the chronic aberration in the baseline levels of TNF in human circulatory system has been attributed to the pathogenesis of numerous diseases, including rheumatoid arthritis, osteoporosis, sepsis and cancer [4, 5].
The vast majority of TNF related biological processes are initiated by the death-domain (DD) containing TNFR1, which is also called TNFRSF1A. Unlike TNFR2, TNFR1 is present in almost all cell types in humans. Upon TNF binding, TNFR1 trimerizes, and its intracellular DD recruits TRADD, which then creates a platform for RIP1 and TRAF2 to collectively form the receptor-signaling complex I. Cellular inhibitor of apoptosis proteins (cIAP)-1 and -2 bind to complex I and, consequently, together with K63-linked ubiquitin chains, modify RIP1 and TRAF2 . This creates docking sites for an E3 ligase or linear ubiquitin chain assembly complex (LUBAC) consisting of heme-oxidized IRP2 ubiquitin ligase-1 (HOIL-1), HOIL-1-interacting protein (HOIP), and SHANK-associated RH domain interacting protein (SHARPIN). Subsequently, the activation of TAK1 and the ubiquitination of NEMO (or IKKγ), a subunit of IKK complex, lead to cell survival or proinflammatory response through NF-κB and MAP kinases activations. Other TRAF superfamily members (TRAF5 and 6) are also known to play a role in the NF-κB and MAP kinases activations [7, 8].
On the other hand, for the apoptotic pathways, clathrin, AP-2 and Dyn first mediate receptor internalization. Receptor-signaling complex I becomes modified, and dissociates from TNFR1, allowing FADD and caspase-8 to form complex II. Within complex II, caspase-8 becomes activated to induce extrinsic apoptosis through caspase-3 activation. Alternatively, caspase-8 activates caspase-7, and eventually, the cleavage of Bid to tBid in the mitochondria activates caspase-9 via cathepsin D. This induces the intrinsic apoptosis through caspase-3 activation.
Due to its ability to signal numerous cellular processes via the survival and death pathways, the TNFR1 signaling research has received immense attention over the years, especially on understanding the downstream signaling cascades to regulate and control proinflammatory diseases and cancer. Despite numerous studies, the control of proinflammatory diseases through therapeutic treatments, where TNF is over-expressed, remains suboptimal. For example, biologic response modifiers or biologics, such as Etanercept and Infliximab, are TNF decoy receptors or antibodies that suppress TNFR1 signaling through competition for TNF. Although these drugs have shown successful downregulation of inflammation in many cases, they can immuno-compromise patients to secondary infections such as tuberculosis , or have been ineffective in a substantial number of administered patients .
To find alternatives, there have been major efforts on selectively suppressing the intracellular signaling of TNFR1. For example, genetic knockouts (KOs) of TRAFs and TRADD acting on the proinflammatory pathways have been investigated [7, 8, 11]. However, the experimental outcomes, so far, have not been optimistic. In TRAF2 KO, there is compensatory activation of NF-κB through TRAF5  or TRAF6 , and vice-versa. On the other hand, TRADD KO almost completely abolishes NF-κB activation , which is not desirable for the general survivability of cells. Thus, a systemic approach where the propagation of signal transduction to all known branching pathways during target intervention should be monitored. This will allow the elucidation of effective target candidate(s) that overcomes and balances the deficiencies of current investigations.
In this paper, we adopted a systems biology approach to study TNFR1 signaling dynamics. Firstly, we developed a computational model of TNF-induced proinflammatory response leading to NF-κB, MAP kinase activations, and three groups of gene expressions (classified according to their temporal profiles ). The model is based on the perturbation-response approach [13–16], which has been successfully used to elucidate novel signaling features and behaviors in Toll-like receptor-4 [17, 18], -3 , and TNF-related apoptosis-inducing ligand (TRAIL) signaling . Secondly, the TNFR1 model parameters were selected to fit the temporal activation profiles of NF-κB and MAP kinase p38 for fibroblast cell type in several available conditions (wildtype , TRAF2 KO , TRAF5 KO , TRAF2/TRAF5 double KO (DKO) , TRAF6 KO , TRADD KO  and RIP1 KO ). Using the resultant TNFR1 model with robust parameters, we performed simulations of multiple in silico KOs to determine an optimal target that suppresses, but not abolishes, proinflammatory genes. Finally, to validate the modeling results, we performed experiments measuring various key proinflammatory gene expressions in MEF and 3T3 cells for TNF stimulation. Overall, our study presents evidence that systems biology research can be useful to elucidate important target(s) to suppress proinflammatory diseases such as rheumatoid arthritis and osteoarthritis.
TNFR1 signaling topology and model
To develop a computational model of proinflammatory TNFR1 signaling dynamics, we first require the known signal transduction pathways. We curated the KEGG database, and performed literature survey of the latest TNF research. After carefully considering several sources, we were able to propose a signaling topology mainly by combining the knowledge from KEGG, Falschlehner et al. (2012) and Wertz et al. (2010) [6, 22] (Figure 1).
Next, to simulate TNF-induced dynamics of NF-κB and MAPK activations using the topology, we developed a dynamic model based on perturbation-response approach (Materials and Methods), using COPASI simulation platform . Unlike common biochemical reaction models [24, 25], the perturbation-response approach does not require detailed knowledge of all signaling species and their reaction kinetics. This is because it analyses the response waves of signal transduction instead of individual reaction kinetics [13–15, 17–20]. The response waves can be approximated using linear response rules (Response Rules, Additional file 1: Figure S1) combined with the law of mass conservation, and this approach has been previously used to successfully model the TLRs and TRAIL signaling pathways [17–20].
Briefly, each reaction in the model is represented by a first-order response equation with activation or deactivation term. The activation term generally refers to protein binding, transformation, complex formation, phosphorylation and transcription. The deactivation term refers to protein unbinding, dephosphorylation and negative regulation such as mRNA decay through microRNA regulation.
Simulating TNF-induced Ν F-κΒ and MAP kinase dynamics
The parameters of the initial model (rate constants, or the elements of Jacobian matrix J, Materials and Methods) were estimated by fitting the simulation profiles with experimental profiles of signaling molecules where data is available. We obtained published semi-quantitative experimental profiles of IκBα phosphorylation (Ν F-κΒ activation) and p38 (MAP kinase) activation in wildtype and various genetically mutant MEFs generally treated with 10 ng/mL of TNF (Figure 2A, Additional file 1: Figure S2 and Table S1) [7, 8, 11, 21]. (Note that the kinetics of other MAP kinases, JNK and ERK, were also similar to p38 [7, 8, 11, 21]. Thus, we used p38 as a representative MAP kinase for our investigation).
The parameter values were selected by using Genetic Algorithm  module in COPASI software  to fit the experimental profiles (Figure 2A, WT). Following, we performed sensitivity analysis (Materials and Methods) of the model parameters and found them to be robust to a small degree of uncertainty to their values (Additional file 1: Table S2). As a further validity of the parameter values, we tested the wildtype model in other conditions, namely TRAF2 KO, TRAF5 KO, TRAF2/5 double KO, TRAF6 KO, RIP1 KO and TRADD KO (Figure 2B). (Note that in silico KOs were generated from the wildtype model by setting the activation parameter value of the KO molecule to null).
Remarkably, we were able to obtain a single set of model parameters (Table 1, reactions 1–29 and see Additional file 2 for the TNFR1 model A in SBML format), which could be used to simulate the semi-quantitative profiles of IκBα phosphorylation and p38 kinase activation in multiple experimental conditions. In wildtype, TRAF2 KO, TRAF5 KO and TRAF6 KO, the IκBα phosphorylation and p38 kinase activation reach peak values around 15 min and gradually decay at 30 min. Notably, TRAF6 KO shows enhanced IκBα phosphorylation and p38 kinase activation due to Signaling Flux Redistribution (Response Rule 5, Additional file 1: Figure S1) . In the remaining conditions, the activation levels of both molecules are very weak (RIP1 KO and TRAF2/5 DKO) or absent (TRADD KO).
It is noteworthy that although there have been previous models on TNF signaling [24, 27, 28], to our knowledge, this is the first time a single model of TNF signaling with fixed parameter values recapitulates the proinflammatory signaling dynamics in multiple experimental conditions.
To compare our linear response model (TNFR1 model A) simulations with other models that contain more detailed descriptions of IKK  and MAPK  signaling, using higher order terms and Michaelis-Menten type kinetics, we developed an alternative TNFR1 model B incorporating the relevant reaction details (Additional file 1: Table S3). Notably, the simulations of TNFR1 models A and B show very similar dynamics for a fixed amount of TNF perturbation (Additional file 1: Figure S3). Thus, we concur that our linear response model can be appropriately used for further investigations.
Simulating distinct TNF-induced gene expression patterns
Next, we extended the TNFR1 model (we will now simply call TNFR1 model A as TNFR1 model) to simulate downstream proinflammatory gene expression dynamics. Recently, time-series high throughput microarray and quantitative real time PCR experiments on TNF simulated mouse 3T3 fibroblasts cells have revealed 3 distinct groups of upregulated gene expression patterns, with possibly corresponding distinct biological roles [12, 30]. The groups were labeled into “early I”, “intermediate or middle II” and “late III” response, according to their time to reach peak expressions between 0.5-1, 2–3, and 6–12 h, respectively, after TNF stimulation (Figure 3A) [12, 30, 31]. Here, we extended the TNFR1 model to simulate the temporal profiles of the 3 groups of gene expressions.
According to our modeling approach, the time to peak activation can be controlled by reaction parameter values and/or the number of signaling intermediates [15, 17–20]. Briefly, decreasing (increasing) the activation or transcription parameter value will show lower (sharper) gradients of formation part of the expression profiles. Alternatively, decreasing (increasing) the deactivation or decay parameter value will show lower (sharper) gradients of depletion part of expression profiles (Response Rule 1, Additional file 1: Figure S1). In addition, inserting intermediary reactions between transcription process and gene induction will increase delay for gene expression dynamics (Response rule 2, Additional file 1: Figure S1). The intermediates can represent the complexities of transcription process involving the pre-initiation, initiation, promoter clearance, elongation and termination , or post-transcriptional processes such as messenger RNA editing and splicing. Using this approach, the TNFR1 model was extended to simulate the temporal dynamics of groups I, II and III genes. Note that the response rules (Additional file 1: Figure S1) are used to modify an initial signaling topology only after all parameter space has been exhaustively searched, and a reasonable model fit is unable to be achieved .
Previous investigations on the 3 groups of genes have indicated distinct mechanisms for the differential dynamical response [12, 30]. Hao and Baltimore have found lesser presence of AU Rich Element (ARE) region on the 3’UTR of group III genes, targeted by microRNAs and ARE-binding proteins (such as tristetraprolin) that enhance RNA decay processes. Hence, it was postulated as one possible reason for the lower decay response of group III genes compared with genes from groups I and II . More recently, by studying the kinetics of pre-mRNA and mRNA, Hao and Baltimore observed delays in splicing of groups II and III genes compared to group I genes. The differential delays were suggested as another biological mechanism for the distinct gene profiles .
In our extended model, we, therefore, considered both mechanisms to reproduce the temporal profiles of the 3 groups of genes. Notably, our simulations of pre-mRNA and mRNA for all groups of genes matched the data of Hao and Baltimore for the first 60 min (Additional file 1: Figure S4). However, subsequently for 12 h, although the simulations of groups I and II genes were recapitulated, group III simulation was poor (Figure 3B, blue line). Specifically, reducing the parameter value for the decay term representing lower miRNA and ARE-binding proteins regulating decay processes (Response rule 1, Additional file 1: Figure S1), and adding intermediates (Response rule 2, Additional file 1: Figure S1) to provide delays in RNA splicing in our model were not sufficient to produce the continuous activation of group III genes (Figure 3C, cyan dotted line).
To overcome the shortfall in the model simulations, we hypothesized that novel activation or transcription term(s) (positive feedback) may be present to provide additional flux for the continuous increase in group III expressions (Response rule 4, Additional file 1: Figure S1). This could result from secondary post-transcriptional/translational mechanisms through i) autocrine signaling such as IL-1 , IL-6  or TGF-β  signaling (Figure 3D), or ii) cytosolic feedback mechanisms specifically for group III genes  (Figure 3E). Thus, a novel feedback mechanism predominantly affecting the transcription of group III genes was added to the TNFR1 model (Table 1, equations 30–65 and Additional file 2).
The modified TNFR1 model with feedback mechanisms to group III genes produced simulations that matched all 3 groups of gene expression profiles (Figure 3A and C, solid lines). To scrutinize the feedback mechanism, we re-monitored the simulation profile of NF-κB for 6 hours (Additional file 1: Figure S5). The resultant profile mimics the damped oscillatory dynamics of NF-κB previously observed in murine fibroblasts . Overall, these data suggest that low miRNA regulation and additional delay in RNA splicing are not sufficient to produce the continuous activation of group III genes, and that a novel transcription process, possibly through secondary post-transcriptional/translational autocrine signaling, such as IL-1 signaling or other novel feedback mechanisms that activate NF-κB, and not MAPK (Additional file 1: Figure S5), are required.
Predicting key target for regulating proinflammatory response
Now that the TNFR1 model is able to successfully simulate the three groups of upregulated genes in wildtype, we investigated the significance and effect of removing or suppressing key intracellular signaling molecules for controlling proinflammatory response, in silico.
It is well known that TNFR1 signaling is enhanced in proinflammatory diseases and cancer [1–4]. To investigate which known molecules would be potential target to regulate the cell survival or proinflammatory activity, we performed in silico KOs of all possible signaling molecules within the TNFR1 model. In total, we simulated groups I, II and III dynamic gene expressions in 12 (TRADD, cIAP1/2, TRAF2, TRAF5, TRAF6, RIP1, SHARPIN, LUBAC, TAK1 complex (TAK1/TAB1/2), IκBα, MKK3/6 and p38) KO conditions and compared with wildtype profiles (Additional file 1: Figure S6).
Among the candidates, the removal of TAK1 complex or RIP1 produced the most noticeable downregulation of all 3 gene groups, which chiefly consist of well-known proinflammatory mediators (Figure 4). However, in TAK1 complex KO, our simulations show almost no induction for group 1 genes. The substantial impairment in gene expressions (> 90%) is usually detrimental to the general survivability of living cells, and this has been particularly demonstrated in TAK1-deficient mice [37, 38]. RIP1, on the other hand, showed about 50-70% impairment compared to wildtype peak expressions. Our simulations, therefore, suggest that RIP1 is possibly a crucial single molecule target for controlling enhanced proinflammatory response due to TNFR1 signaling in proinflammatory disease conditions, such as in rheumatoid arthritis, without compromising the normal functioning of other cellular activities.
Experimental inhibition of RIP1 downregulates proinflammatory genes in TNF stimulation
To verify the predictions of TNFR1 model simulations, we prepared corresponding MEF and BALB/3T3 cells treated with TNF in wildtype and in RIP1 suppression. Necrostatin-1 (Nec-1) was originally identified as a potent small molecule inhibitor of necroptosis or non-apoptotic cell death . Further interests in Nec-1 led to its specificity towards the inhibition of RIP1 . Although Nec-1 has recently been extensively studied, its effect on the expressions of groups I, II and III genes in TNF stimulation remains largely unknown. Therefore, here, we used Nec-1 to suppress RIP1 in vivo.
To check the effect of cell death by Nec-1, we compared MEF and BALB/3T3 cells treated with different doses of Nec-1 in the presence or absence of TNF (Additional file 1: Figure S7). The data revealed that Nec-1 has no substantial effect on cell death after 24 h incubation, and hence, could be tested for its efficacy on the 3 groups of genes. We next performed quantitative RT-PCR for a total of 10 genes: Il6, Tnfaip3, Jun, Nfkbia (group I), Ccl7, Vcam1, Cxcl10 (group II), and Mmp3, Mmp13, Enpp2 (group III). We intentionally included key proinflammatory mediators, genes of matrix metalloproteinase (Mmp3, Mmp13), which are known to degrade collagen in cartilage and thereby enhance rheumatoid arthritis and osteoarthritis progression [41–44].
A previous study has shown that 30 μM of Nec-1 effectively inhibited RIP1 kinase activity . Therefore, we investigated gene expressions for cells stimulated with 10 ng/mL TNF, in the presence or absence of 30 μM Nec-1 for a period of 10 hours with measurements made at least every hour (Figure 5). Remarkably, as predicted by the TNFR1 model, RIP1 inhibition by Nec-1 resulted in the suppression of all 3 groups of genes. The effect of suppressing RIP1 is significant for groups I and II genes in both MEF and BALB/3T3 cells, especially during the first 2–3 hours after stimulation. For group III genes, Nec-1 had more pronounced effect in MEF compared with BALB/3T3 cells. Overall, these results are consistent with the TNFR1 model predictions that suppressing RIP1 in TNF stimulation significantly impairs the activation of all 3 groups of genes.
TNF is a crucial cytokine that regulates myriad vital cellular processes. However, its levels are enhanced in major proinflammatory diseases. Here, to understand the TNF-induced proinflammatory signaling process, and to carefully regulate its dynamic response, a systems biology approach was adopted. We first developed a dynamic computational model using well-established publicly available experimental data of NF-κB, MAP kinase p38, and the average profiles of 3 groups of 180 upregulated genes in mouse fibroblast cells.
Despite the simplicity of using first-order response equations to simulate the profiles of the intracellular molecules, the computational model of TNFR1 recapitulated the experimental response in wildtype and several mutant conditions for NF-κB and p38 activations. This result is surprising, as we know that the innate immune response of TNF is highly complex. It is important to note here that there have been previous other computational efforts on NF-κB and MAPK signaling that had utilized detailed biochemical reactions modeling, to elucidate local properties of signal transduction, such as the ability of common molecules to produce distinct feedback mechanisms to different stimuli [24, 45, 46]. In our work, however, we have shown that even a simpler representation of the signal transduction pathways, through first order response equations and the law of mass conservation can reproduce experimental dynamics. This strongly indicates the presence of simple organizing rules governing the deterministic population average signaling response [47–52].
Next, through the analyses of downstream temporal gene expression profiles, the model suggests the presence of additional novel post-transcriptional/translational processes that is required for the continuous activation of group III genes. This result is additional to previous postulations, which had indicated that the continuous activation is due to lesser ARE region for group III genes leading to a very low decay process , and due to the presence of differential delays in the RNA splicing process . Our model suggests that, on top of these effects, a novel time-delayed secondary transcriptional mechanism is required.
Literature survey indicates that the novel positive feedback processes could be a result of autocrine signaling, example through IL-1 or IL-6, or derive from a still unknown intracellular feedback mechanisms regulating mainly the promoter regions of group III genes. For example, the role of interferon regulatory factor (IRF) family in inducing Ccl5 or RANTES expression, which belongs to one of the group 3 genes, is reported in a previous study , however, was not considered in the initial TNFR1 model. It is, therefore, necessary to perform further experimental work to confirm and elucidate the exact mechanisms for the continuous activations of group III genes.
On the other hand, for down-regulating TNF signaling, which is enhanced in several proinflammatory diseases and cancer, we performed the simulations for 12 in silico KOs of signaling molecules. The resultant simulations indicated that RIP1 is a major regulator of the 3 groups of upregulated gene expressions. To verify the result, we performed experiments on MEF and BALB/3T3 cells using Nec-1 as an inhibitor of RIP1. The measurement of 10 genes belonging to groups I (Il6, Tnfaip3, Jun, Nfkbia), II (Ccl7, Vcam1, Cxcl10) and III (Mmp3, Mmp13, Enpp2) all showed significant impairment with Nec-1 compared to wildtype.
Most importantly, the expressions of key proinflammatory genes such as Il6, Vcam1, Ccl7, Mmp3, Mmp13, enhanced in rheumatoid arthritis and osteoarthritis [41, 44], were reduced. In particular are the levels of matrix metalloproteinase genes Mmp3, Mmp13, which are known to directly affect type II collagen in bone cartilages and degrade the extracellular matrix. Although recent therapeutics have been focusing on the specific regulations of MMPs [42–44, 54], it remains to be seen what effect such treatments will have on other proinflammatory or vital genes.
In summary, our approach provides a systemic analysis of TNFR1 signaling, and suggests Nec-1 is potentially an important therapeutic target for effectively regulating major proinflammatory mediators in chronic diseases where TNF is overexpressed.
Materials and methods
The model is based on perturbation-response approach [15, 17–20]. The basic principle behind the approach is to induce a controlled perturbation of input reaction species of a system (TNFR1), and monitor the response of the activation/concentration levels of other output species (e.g. TAK1, p38, NF-κB, Il6, etc.) from steady-state. To briefly explain the principle, let a stable network consisting of n species be perturbed from the reference steady-state. In general, the resultant changes in the concentration of species are governed by the kinetic evolution equation [13, 14]:
where the corresponding vector form of equation 1 is . F is a vector of any non-linear function including diffusion and reaction of the species vector X = (X 1, X 2 , .., X n ), which represents activated concentration levels of reaction species. The response to perturbation can be written by X = X0 + δX, where X0 is the reference steady-state vector and δX is the relative response from steady-states (δX t=0 = 0).
The generally non-linear kinetic evolution equation 1 can be approximated or linearized by using Taylor series:
As the general volume of perturbing substance is usually very small (order of 1%) compared to the total volume of cells that are perturbed , now consider a small perturbation around the steady-state in equation 2, in which higher-order terms become negligible and result in the approximation of the first-order term. In vector form (note the change from partial derivative to total derivative of time), where the zeroth order term F(X0) = 0 at the steady-state X0 and the Jacobian matrix, or linear stability matrix, is . The elements of J, based on the initial activation topology, are chosen by fitting δX with corresponding experimental profiles. Hence, the amount of response (flux propagated) along a biological pathway can be approximated using first order mass-action response, i.e. . That is, the basic principle so far suggests that the response rate of species in a mass-conserved system at an initial steady-state can be approximated by first order mass-action response equations, given a small perturbation to one or more species.
Note that Jacobian matrix elements (or response coefficients) can include not only reaction information, but also spatial information such as diffusion and transport mechanisms. Thus, each species in the perturbation-response model can represent a molecule, a different modified state of a molecule (e.g. ubiquitinated state) or a molecular process such as diffusion, endocytosis, etc. That is, each species in the biological network does not necessarily represent a specific molecular species. For illustration, in a pathway X 1 → X 2 → X 3 → X 4 → X 5 , X 1 to X 5 can each be a different species or the same species at different stages in signaling, for example, X 1 being internalized (becoming X 2 ), transported to a different organelle (X 3 ), ubiquitinated (X 4 ) and become part of a protein complex (X 5 ).
The complete SBML version of TNFR1 Models A & B are available in Additional file 2.
We performed a sensitivity analysis to test the robustness of the optimized model parameters using the COPASI sensitivities module with default values. The variation in the response of signaling molecules/steps, x i (t), was analyzed when a small variation of each model parameter k j was applied. The response sensitivity coefficient  of the ith molecule with regard to the jth parameter is defined by
The obtained values, R i,j are then scaled, to reflect the relative changes in response, such as a change of p% in the value of parameter k j , results in a R i , j ·p% change in the value of the peak activation of the ith molecule. The response sensitivity coefficients of p38, IκBα, and groups I, II and III genes were obtained at peak time (t = 15 min for p38 and IκB, 30 min, 2 h and 12 h for groups I, II and III respectively, see Additional file 1: Table S2).
Reagents and cell culture
Recombinant mouse TNF was purchased from R&D systems. Necrostatin-1 was purchased from Merck Millipore. 3T3 cells were obtained from JCRB cell bank. 3T3 and MEF were grown in DMEM (Nissui Seiyaku Co.) containing 10% calf serum, 100 U/mL of penicillin at 37°C in a 5% CO2 humidified atmosphere.
Evaluation of cell survival by 3-(4,5-dimethylthiazol-2-yl)- 2,5-diphenyltetrazolium bromide (MTT) assay
The sensitivity of cells to hyperosmotic stress was measured with the MTT colorimetric assay in 96-well plates. Cells (2 × 104) were inoculated in each well and incubated for 24 h. Thereafter, 50 μL of MTT (2 mg/mL in PBS) was added to each well and the plates were incubated for a further 2 h. The resultant formazan was dissolved with 100 μL of dimethyl sulfoxide after aspiration of culture medium. Plates were placed on a plate shaker for 1 min and then read immediately at 570 nm using TECAN microplate reader with Magellan software (Männedorf, Switzerland).
Real-time PCR analysis
Total cellular RNA was extracted from cells using the TRIzol reagent according to the manufacturer’s instructions (Invitrogen). One microgram of RNA was reverse-transcribed using a first-strand cDNA synthesis kit (ReverTra Aceα; Toyobo). Quantitative real-time PCR was performed using SYBR premix Ex Taq (Takara) on the Applied Biosystems StepOnePlusTM according to the technical brochure of the company. RT-PCR primers designed in this study are listed in Additional file 1: Table S4. Quantitative measurements were determined using the ΔΔCt method and expression of GAPDH was used as the internal control. Melt curve analyses of all real-time PCR products were performed and shown to produce the sole DNA duplex.
Tumor necrosis factor
Murine embryonic fibroblast
Receptor-interacting protein 1
Tumor necrosis factor receptor 1 associated death domain protein
Cellular inhibitor of apoptosis proteins
Linear ubiquitin chain assembly complex
Heme-oxidized iron regulatory protein 2 ubiquitin ligase-1
SH3 and multiple ankyrin repeat domains protein-associated RH domain interacting protein
Transforming growth factor β(TGFβ)-activated kinase 1
fas-associated death domain protein
TNF-related apoptosis-inducing ligand
Double knock out
c-Jun N-terminal kinases
Extracellular signal-regulated kinase
Inhibitors of NF-κB
AU rich element
Reverse transcription polymerase chain reaction
Interferon regulatory factor
O’Malley WE, Achinstein B, Shear MJ: Action of bacterial polysaccharide on tumours: II: damage of sarcoma 37 by serum of mice treated with serratia marcescens polysaccharide, and induced tolerance. J Natl Cancer Inst. 1962, 29: 1169-1175.
Balkwill F: Tumour necrosis factor and cancer. Nat Rev Cancer. 2009, 9: 361-371. 10.1038/nrc2628.
Locksley RM, Killeen N, Lenardo MJ: The TNF and TNF receptor superfamilies: integrating mammalian biology. Cell. 2001, 104: 487-501. 10.1016/S0092-8674(01)00237-9.
Bradley JR: TNF-mediated inflammatory disease. J Pathol. 2008, 214: 149-160. 10.1002/path.2287.
van Horssen R, Ten Hagen TL, Eggermont AM: TNF-alpha in cancer treatment: molecular insights, antitumor effects, and clinical utility. Oncologist. 2006, 4: 397-408.
Falschlehner C, Boutros M: Innate immunity: regulation of caspases by IAP-dependent ubiquitylation. EMBO J. 2012, 31: 2750-2752. 10.1038/emboj.2012.148.
Tada K, Okazaki T, Sakon S, Kobarai T, Kurosawa K, Yamaoka S, Hashimoto H, Mak TW, Yagita H, Okumura K, Yeh WC, Nakano H: Critical roles of TRAF2 and TRAF5 in tumor necrosis factor-induced NF-κB activation protection from cell death. J Biol Chem. 2001, 276: 36530-36534. 10.1074/jbc.M104837200.
Funakoshi-Tago M, Kamada N, Shimizu T, Hashiguchi Y, Tago K, Sonoda Y, Kasahara T: TRAF6 negatively regulates TNFalpha-induced NF-kappaB activation. Cytokine. 2009, 45: 72-79. 10.1016/j.cyto.2008.10.010.
Fallahi-Sichani M, Flynn JL, Linderman JJ, Kirschner DE: Differential risk of tuberculosis reactivation among anti-TNF therapies is due to drug binding kinetics and permeability. J Immunol. 2012, 188: 3169-3178. 10.4049/jimmunol.1103298.
Wiens A, Venson R, Correr CJ, Otuki MF, Pontarolo R: Meta-analysis of the efficacy and safety of Adalimumab, Etanercept, and Infliximab for the treatment of rheumatoid arthritis. Pharmacotherapy. 2010, 30: 339-353. 10.1592/phco.30.4.339.
Ermolaeva MA, Michallet MC, Papadopoulou N, Utermöhlen O, Kranidioti K, Kollias G, Tschopp J, Pasparakis M: Function of TRADD in tumor necrosis factor receptor 1 signaling and in TRIF-dependent inflammatory responses. Nat Immunol. 2008, 9: 1037-1046. 10.1038/ni.1638.
Hao S, Baltimore D: The stability of mRNA influences the temporal order of the induction of genes encoding inflammatory molecules. Nat Immunol. 2009, 10: 281-288.
Ross J: Determination of complex reaction mechanisms: analysis of chemical, biological and genetic networks. J Phys Chem A. 2008, 112: 2134-2143. 10.1021/jp711313e.
Moran F, Vlad MO, Bustos M, Trivino JC, Ross J: Species connectivities and reaction mechanisms from neutral response experiments. J Phys Chem A. 2007, 111: 1844-1851. 10.1021/jp0661793.
Selvarajoo K, Tomita M, Tsuchiya M: Can complex cellular processes be governed by simple linear rules?. J Bioinform Comput Biol. 2009, 7: 243-268. 10.1142/S0219720009003947.
Vance W, Arkin A, Ross J: Determination of causal connectivities of species in reaction networks. Proc Natl Acad Sci USA. 2002, 99: 5816-5821. 10.1073/pnas.022049699.
Selvarajoo K: Discovering differential activation machinery of the toll-like receptor 4 signaling pathways in MyD88 knockouts. FEBS Lett. 2006, 580: 1457-1464. 10.1016/j.febslet.2006.01.046.
Selvarajoo K, Takada Y, Gohda J, Helmy M, Akira S, Tomita M, Tsuchiya M, Inoue J, Matsuo K: Signaling flux redistribution at toll-like receptor pathway junctions. PLoS ONE. 2008, 3: e3430-10.1371/journal.pone.0003430. DOI:10.1371/journal.pone.0003430
Helmy M, Gohda J, Inoue J, Tomita M, Tsuchiya M, Selvarajoo K: Predicting novel features of Toll-like receptor 3 signaling in macrophages. PLoS ONE. 2009, 4: e4661-10.1371/journal.pone.0004661. DOI: 10.1371/journal.pone.0004661
Piras V, Hayashi K, Tomita M, Selvarajoo K: Enhancing apoptosis in TRAIL-resistant cancer cells using fundamental response rules. Sci Rep. 2011, 1: 144-DOI:10.1038/srep00144
Devin A, Cook A, Lin Y, Rodriguez Y, Kelliher M, Liu Z: The distinct roles of TRAF2 and RIP in IKK activation by TNF-R1: TRAF2 recruits IKK to TNF-R1 while RIP mediates IKK activation. Immunity. 2000, 12: 419-429. 10.1016/S1074-7613(00)80194-6.
Wertz IE, Dixit VM: Regulation of death receptor signaling by the ubiquitin system. Cell Death Differ. 2010, 17: 14-24. 10.1038/cdd.2009.168.
Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI–a COmplex PAthway SImulator. Bioinformatics. 2006, 22: 3067-3074. 10.1093/bioinformatics/btl485.
Werner SL, Barken D, Hoffmann A: Stimulus specificity of gene expression programs determined by temporal control of IKK activity. Science. 2005, 309: 1857-1861. 10.1126/science.1113319.
Tasseff R, Nayak S, Salim S, Kaushik P, Rizvi N, Varner JD: Analysis of the molecular networks in androgen dependent and independent prostate cancer revealed fragile and robust subsystems. PLoS One. 2010, 5: e8864-10.1371/journal.pone.0008864. DOI:10.1371/journal.pone.0008864
Carroll DA: Chemical laser modelling with genetic algorithms. AIAA. 1996, 34: 338-346. 10.2514/3.13069.
Werner SL, Kearns JD, Zadorozhnaya V, Lynch C, O’Dea E, Boldin MP, Ma A, Baltimore D, Hoffmann A: Encoding NF-κB temporal control in response to TNF: distinct roles for the negative regulators IκBα and A20. Genes Dev. 2008, 22: 2093-2101. 10.1101/gad.1680708.
Cho KH, Shin SY, Lee HW, Wolkenhauer O: Investigations into the analysis and modeling of the TNF alpha-mediated NF-kappa B-signaling pathway. Genome Res. 2003, 13: 2413-2422. 10.1101/gr.1195703.
Kholodenko BN: Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem. 2000, 267: 1583-1588. 10.1046/j.1432-1327.2000.01197.x.
Hao S, Baltimore D: RNA splicing regulates the temporal order of TNF-induced gene expression. Proc Natl Acad Sci U S A. 2013, 110: 11934-11939. 10.1073/pnas.1309990110.
Tian B, Nowak DE, Brasier AR: A TNF-induced gene expression program under oscillatory NF-kappaB control. BMC Genomics. 2005, 6: 137-154. 10.1186/1471-2164-6-137.
Solomon EP, Berg LR, Martin DW: Biology. 2011, Thomson Brooks/Cole: Pacific Grove, 9
Chaudhry SI, Hooper S, Nye E, Williamson P, Harrington K, Sahai E: Autocrine IL-1beta-TRAF6 signalling promotes squamous cell carcinoma invasion through paracrine TNFalpha signalling to carcinoma-associated fibroblasts. Oncogene. 2013, 32: 747-758. 10.1038/onc.2012.91.
Grivennikov S, Karin M: Autocrine IL-6 signaling: a key event in tumorigenesis?. Cancer Cell. 2008, 13: 7-9. 10.1016/j.ccr.2007.12.020.
Ihn H: Autocrine TGF-beta signaling in the pathogenesis of systemic sclerosis. J Dermatol Sci. 2008, 49: 103-113. 10.1016/j.jdermsci.2007.05.014.
Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IκB–NF-κB signaling module: temporal control and selective gene activation. Science. 2002, 298: 1241-1245. 10.1126/science.1071914.
Tang M, Wei X, Guo Y, Breslin P, Zhang S, Zhang S, Wei W, Xia Z, Diaz M, Akira S, Zhang J: TAK1 is required for the survival of hematopoietic cells and hepatocytes in mice. J Exp Med. 2008, 205: 1611-1619. 10.1084/jem.20080297.
Lamothe B, Lai Y, Xie M, Schneider MD, Darnay BG: TAK1 is essential for osteoclast differentiation and is an important modulator of cell death by apoptosis and necroptosis. Mol Cell Biol. 2013, 33: 582-595. 10.1128/MCB.01225-12.
Degterev A, Huang Z, Boyce M: Chemical inhibitor of nonapoptotic cell death with therapeutic potential for ischemic brain injury. Nat Chem Biol. 2005, 1: 112-119. 10.1038/nchembio711.
Degterev A, Hitomi J, Germscheid M, Ch’en IL, Korkina O, Teng X, Abbott D, Cuny GD, Yuan C, Wagner G, Hedrick SM, Gerber SA, Lugovskoy A, Yuan J: Identification of RIP1 kinase as a specific cellular target of necrostatins. Nat Chem Biol. 2008, 4: 313-321. 10.1038/nchembio.83.
Liacini A, Sylvester J, Li WQ, Zafarulla M: Inhibition of interleukin-1-stimulated MAP kinases, activating protein-1 (AP-1) and nuclear factor kappa B (NF-kappa B) transcription factors down-regulates matrix metalloproteinase gene expression in articular chondrocytes. Matrix Biol. 2002, 21: 251-262. 10.1016/S0945-053X(02)00007-0.
Liacini A, Sylvester J, Li WQ, Huang W, Dehnade F, Ahmad M, Zafarullah M: Induction of matrix metalloproteinase-13 gene expression by TNF-alpha is mediated by MAP kinases, AP-1, and NF-kappaB transcription factors in articular chondrocytes. Exp Cell Res. 2003, 288: 208-217. 10.1016/S0014-4827(03)00180-0.
Sellam J, Berenbaum F: The role of synovitis in pathophysiology and clinical symptoms of osteoarthritis. Nat Rev Rheumatol. 2010, 6: 625-635. 10.1038/nrrheum.2010.159.
Roman-Blas JA, Jimenez MD: NF-κB as a potential therapeutic target in osteoarthritis. Osteoarthritis Cartilage. 2006, 14: 839-848. 10.1016/j.joca.2006.04.008.
Santos SD, Verveer PJ, Bastiaens PI: Growth factor–induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat Cell Biol. 2007, 9: 324-330. 10.1038/ncb1543.
Bruggeman FJ, Westerhoff HV, Hoek JB, Kholodenko BN: Modular response analysis of cellular regulatory networks. J Theor Biol. 2002, 218: 507-520.
Selvarajoo K, Tomita M: Physical laws shape biology. Science. 2013, 339: 646-
Selvarajoo K: Interpreting the dynamics and patterns of living systems. Bioscience. 2013, 63: 721-722.
Selvarajoo K: Macroscopic law of conservation revealed in the population dynamics of Toll- like receptor signaling. Cell Commun Signal. 2011, 9: 9-10.1186/1478-811X-9-9.
Selvarajoo K: Understanding multimodal biological decisions from single cell and population dynamics. WIREs Syst Biol Med. 2012, 4: 385-399. 10.1002/wsbm.1175.
Selvarajoo K: Uncertainty and certainty in cellular dynamics. Front Genet. 2013, 4: 68-DOI:10.3389/fgene.2013.00068
Selvarajoo K, Giuliani A: Finding self-organization from the dynamic gene expressions of innate immune responses. Front Physiol. 2012, 3: 192-DOI:10.3389/fphys.2012.00192
Yarilina A, Park-Min KH, Antoniv T, Hu X, Ivashkiv LB: TNF activates an IRF1- dependent autocrine loop leading to sustained expression of chemokines and STAT1-dependent type I interferon-response genes. Nat Immunol. 2008, 9: 378-387.
Kaneva MK, Kerrigan MJ, Grieco P, Curley GP, Locke IC, Getting SJ: Chondroprotective and anti-inflammatory role of melanocortin peptides in TNF-alpha activated human C-20/A4 chondrocytes. Br J Pharmacol. 2012, 167: 167-179.
Theobald U, Mailinger W, Baltes M, Rizzi M, Reuss M: In vivo analysis of metabolic dynamics in Saccharomyces cerevisiae: I: Experimental observations. Biotechnol Bioeng. 1997, 55: 305-316. 10.1002/(SICI)1097-0290(19970720)55:2<305::AID-BIT8>3.0.CO;2-M.
Zi H: Sensitivity analysis approaches applied to systems biology models. IET Syst Biol. 2011, 5: 336-346. 10.1049/iet-syb.2011.0015.
We thank Kiyotoshi Sato, Tomoyoshi Soga and Mitsuhiro Kitagawa for experimental support.
This work was supported by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research J13108 (K.S.), Scientific Research F13804 (K.H.), and Tsuruoka City, Yamagata Prefecture.
The authors declare that they have no competing interests.
KH and KS conceptualized and designed the study. KH and ST performed the wet experiments. VP and KS performed the computational simulations. MT provided cells, reagents and discussions. KH and KS wrote the paper. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1: Figure S1: Response rules. Figure S2. Experimental raw data used for model fitting. Figure S3. Experimental vs. simulated profiles of IκBα and p38 activations in wildtype and mutant conditions using TNFR1 model B. Figure S4. Simulation of pre-mRNA and mRNA expression profiles of groups I, II and III genes. Figure S5. Simulation of NF-κB activation profiles with and without feedback mechanisms. Figure S6. The effects of in silico KOs on the expression profiles of groups I, II and III genes. Table S1. Estimation of the relative intensities of IκBα and p38 activation dynamics. Table S2. Sensitivity analysis of TNFR1 model A. Table S3. TNFR1 model B details. Table S4. List of primer sequences for RT-PCR. (DOCX 7 MB)
About this article
Cite this article
Hayashi, K., Piras, V., Tabata, S. et al. A systems biology approach to suppress TNF-induced proinflammatory gene expressions. Cell Commun Signal 11, 84 (2013). https://doi.org/10.1186/1478-811X-11-84
- Cell signaling
- Computational model