Predicting Novel Features of Toll-Like Receptor 3 Signaling in Macrophages

The Toll-like receptor (TLR) 3 plays a critical role in mammalian innate immune response against viral attacks by recognizing double-stranded RNA (dsRNA) or its synthetic analog polyinosinic-polycytidylic acid (poly (I∶C)). This leads to the activation of MAP kinases and NF-κB which results in the induction of type I interferons and proinflammatory cytokines to combat the viral infection. To understand the complex interplay of the various intracellular signaling molecules in the regulation of NF-κB and MAP kinases, we developed a computational TLR3 model based upon perturbation-response approach. We curated literature and databases to determine the TLR3 signaling topology specifically for murine macrophages. For initial model creation, we used wildtype temporal activation profiles of MAP kinases and NF-κB and, for model testing, used TRAF6 KO and TRADD KO data. From dynamic simulations we predict i) the existence of missing intermediary steps between extracellular poly (I∶C) stimulation and intracellular TLR3 binding, and ii) the presence of a novel pathway which is essential for JNK and p38, but not NF-κB, activation. Our work shows activation dynamics of signaling molecules can be used in conjunction with perturbation-response models to decipher novel signaling features of complicated immune pathways.


Introduction
Toll-Like Receptors (TLRs) play a major role in innate immunity, the first line of mammalian defense against invading pathogens and crucial for antigen-specific acquired immunity development [1,2]. Upon the recognition of pathogen-associated molecular patterns (PAMPs), such as bacterial lipopolysaccharide (LPS) and viral dsRNA, the 13 currently known TLRs trigger predominantly the activation of MAP kinases and several key transcription factors including nuclear factor-kB (NF-kB), activator protein (AP)-1 and interferon regulatory factor (IRF)-3 and 7. This results in the induction of numerous proinflammatory cytokines and type I interferons [3,4,5]. The dysregulation of TLR signaling, therefore, has been attributed to the pathogenesis of major pro-inflammatory illnesses such as the autoimmune diseases [6,7].
TLR3, one of the 4 known intracellular members of TLR family, recognizes dsRNA and poly (I:C) [3,5,8] triggers an innate response independent of the adaptor protein Myeloid Differentiation factor 88 (MyD88), which is required for all other TLRs [8,9]. The specificity of TLR3 response is possibly due to the occurrence of an alanine residue in a critical region of its cytoplasmic domain unlike the proline residue utilized by MyD88 found in other TLRs [5]. Thus, TLR3 initiates its response depending only on the adaptor protein TIR domain-containing adapter-including interferon-b (TRIF) [2,10]. The recruitment of TRIF mediates the signaling process through the activation of key transcription factors NF-kB, AP-1, IRF-3 and 7 [3,5,8]. Although the signaling molecules and their cascades have been broadly investigated, the dynamic outcome of signal transduction between wildtype and genetic mutations still remains poorly understood.
Here, we began the investigation of TLR3 pathway by literature/database curation of the signaling topology in murine macrophages. Next, we analyzed the temporal experimental data of wildtype, TNF Receptor Associated Factor (TRAF)-6 and NFRSF1A-associated via death domain (TRADD)-deficient murine macrophages with poly (I:C) stimulation [11][12] by developing a computational TLR3 model based on perturbation-response approach. This approach does not require the detailed reaction kinetics of each reaction in the signaling topology (as in bottom-up approaches), but rather, considers the activation of signaling molecules as linear response events [13]. Similar modeling approaches have been previously used to infer important biological network features; inferring feedback control of IKK activity in tumour-necrosis factor (TNF) stimulation [14], uncovering switching behaviour of MAPK signaling between epidermal growth factor (EGF) and neuronal growth factor (NGF) stimuli [15], detecting connectivities of reaction molecules [16], predicting missing molecules in TLR4 signaling [17] and signaling flux redistribution (SFR) at pathway junctions [18].
Our TLR3 model simulations were compared with temporal experimental data of NF-kB, JNK and p38 in three experimental conditions; wildtype, TRAF6-deficient and TRADD-deficient macrophages. Collectively, the results suggest i) the existence of novel intermediary steps (e.g. missing cellular processes, proteins or phosphorylation states) between extracellular poly (I:C) stimulation and intracellular TLR3 binding, and ii) the presence of a novel pathway which is essential for JNK and p38 activation.

Results and Discussion
Determination of the TLR3 signaling topology in macrophages TLR3 is expressed in several cell types including macrophages, murine embryonic fibroblasts (MEFs) and dendritic cells (DCs), however, it was not found in B Cells, T Cells and NK cells [8,19,20,21,22]. Moreover, the experimental observations of TLR3 signaling under various genetic knock-outs show controversial roles of certain molecules in different cell types. For instance, Gohda et al. [11] showed that TRAF6 is dispensable for TLR3 induced NF-kB in poly (I:C) stimulated macrophages, while Jiang et al. [23] showed the requirement of TRAF6 to NF-kB activation in MEFs. Also, for NF-kB and MAP kinases activation, TRADD is not critical in macrophages, whereas, it is important for MEFs [12,24]. These data indicate that there is no common topology for TLR3 signaling and, therefore, we cannot combine data obtained from different cell types to create unified TLR3 signaling topology. Instead, independent cell type analysis should be performed.
Here, we investigate the signaling topology for murine macrophages only. It is known that TRIF interacts with TLR3 at the TIR domain and TRAF6 binds to the N-terminal of TRIF [8,11,25]. However, we found, in murine macrophages, the role of TRAF6 in TLR3 signaling is dispensable [11]; the temporal experimental profiles of MAP kinases (JNK and p38) and NF-kB activation to poly (I:C) stimulation in TRAF6 KO were only slightly reduced compared with wildtype levels ( Figure 1A-C WT and TRAF6 KO).
The prediction of missing intermediary cellular processes in poly (I:C) stimulated macrophages So far, we have built the signaling topology for TLR3 pathways in macrophages. Next, to understand the complex dynamic interplay of the various intracellular signaling molecules in the regulation of NF-kB, JNK and p38 in poly (I:C) stimulation, we developed a computational model of the TLR3 signaling (see Materials and Methods). Each signaling response (i.e., the rate of activation of each signaling molecule) in the model is represented by ddX dt~J dX, where dX is relative activated concentration of signaling molecules and the parameters (elements of J) are chosen to fit the semi-quantitative experimental profiles (e.g, western blots, EMSA, etc.) of NF-kB, JNK and p38 of wildtype macrophages stimulated with poly (I:C) ( Figure 1A-C WT). As mentioned earlier, TRAF6 KO does not noticeably affect NF-kB, JNK or p38 activation compared to wildtype ( Figure 1A-C).
Hence, in the model, we used a low parameter value between TRIF and TRAF6 to limit the response flux propagation through TRAF6 such that the removal of TRAF6 in the model will only slightly affect its downstream reactions (see Materials and Methods for details).
We, next, simulated activation of NF-kB, JNK and p38 in wildtype. Our results show the time to reach peak values were 15-  [11,12]. The activation levels were quantified from the western blots using ImageJ [http://rsbweb.nih.gov/ij/]. The x-axis represents the time in minutes and the y-axis represents the relative activation profile. Note: As TRADD KO data is unavailable at 10 min (earliest at 15 min), we could not observe delayed activation of as noted for WT and TRAF6 KO. Therefore, we used dotted line to connect 0 min and 15 min time points. doi:10.1371/journal.pone.0004661.g001 20 min earlier than in actual experiments ( Figure 2B-D). Clearly, this model could not explain the delayed experimental activation of NF-kB, JNK and p38. We can consider time delay processes of signaling events are due to missing molecules/complex formation or spatial movement of molecules [13]. The deterministic kinetic evolution equation used in our model, non-linear F in general (Eq.1, Materials and Methods), can include such information by setting total derivative of time to partial derivative in time and space. In other words, Jacobian matrix J can contain temporal and spatial information of the network process. Since we are not performing spatial simulation, time delay response can be lumped as missing molecules/processes in the network. For example, using a TLR4 model, we previously predicted the delayed activation of TRIF-dependent pathways by using a number of additional response reactions representing missing signaling features (molecules/processes, spatial movements or complex formation) in the original network [17]. Our result was later substantiated by McGettick et al., 2006 [28] who demonstrated two novel signaling molecules (PKCe and TRAM) act upstream of TRIF and, recently, by Kagan et al., 2008 [29] who discovered that the internalization of TLR4 into the endosome is required prior to TRIF-dependent pathway activation. Thus, the delayed activation of NF-kB and MAP Kinases in poly (I:C) stimulation could also be due to missing molecules or cellular processes.
To investigate and locate the likely position of the necessary missing intermediary steps in the TLR3 pathway predicted by the Figure 2. Analysis of TLR3 pathway in macrophages. A) Schematic representation of the determined TLR3 pathway topology in macrophages. dsRNA or poly (I:C) stimulated TLR3 triggers TRIF dependant response by the recruitment of TRIF to the cytoplasmic domain of the receptor which then allows RIP1, TRAF6, TBK1 and TRAF3 to bind with TRIF. This results in the activation of MAP kinases (MKK1/2, MKK3/6 and MKK4/7) and IkB kinase complex; MKK1/2, MKK3/6 and MKK4/7 activate ERK, JNK and p38, respectively and IkBa degradation releases NF-kB. TBK1 phosphorylates IRF-3 and 7. ERK, JNK and p38 translocate to the nucleus and activate the transcription factor AP-1, and NF-kB, IRF-3 and IRF-7 translocate to the nucleus. AP-1 and NF-kB bind to the promoter regions of cytokine genes such as Tnf and Il6 while IRF-3, IRF-7 together with NF-kB bind to the promoter region of chemokine genes such as Cxcl10 and Ccl5 and induce their transcription. Protein-protein interactions between molecules at the two signaling branches analyzed are highlighted in brown and blue. The dotted lines indicate weak activation (see maintext). B), C) and D) show simulations of NF-kB, JNK and p38 activation, respectively, in wildtype (WT). The x-axis represents the time in minutes and the y-axis represents the relative activation profile. doi:10.1371/journal.pone.0004661.g002 model, we further surveyed the literature. TRIF directly bind to the TIR domain of TLR3 and does not require TRAM for its activation [30]. Thus, it is unlikely that missing intermediary steps exist between TLR3 and TRIF. Furthermore, TLR3 KO and TRIF KO both showed similar response; abolished activation of NF-kB and MAP Kinases [30]. These data led us to hypothesize the missing intermediary steps (e.g. signaling molecules or processes) are upstream of TLR3 and we updated our model to begin simulation not from TLR3, but rather from poly (I:C) downwards ( Figure 3A). We also represented each new uncharacterized molecule/process as a signaling intermediate in our model.
To obtain the delayed activation of NF-kB and MAP kinases in accordance with experimental data of wildtype, that is, null activation till 10 min, peak values around 30 min and reduced activation after 60 min for all molecules ( Figure 1A-C), we were required to add three intermediary steps (signaling intermediates) ( Figure 3B-D WT). Having less or more intermediary steps resulted in peak value reaching faster or slower than experimental peak, respectively (data not shown).
While preparing this manuscript, Liu et al. demonstrated that TLR3-ectodomains dimerizes before signal propagation [31]. Our prediction of novel intermediary steps in upstream of TLR3 activation indicates that TLR3-ectodomains dimerization could be one of the missing intermediary steps. The other intermediary steps could indicate the endocytosis of poly (I:C) and its subsequent transport mechanisms to be recognized by TLR3.
The existence of missing pathway for MAP Kinases activation in poly (I:C) stimulation To further test the predictive capability of our model, we simulated TRAF6 KO and TRADD KO (see Materials and Methods). We wondered whether the same TLR3 model could successfully simulate the activation of NF-kB, JNK and p38 in all three conditions: wildtype, TRAF6 KO and TRADD KO. To create KO condition from wildtype model, we set the reactions of the KO molecule null, while all other model parameters retain their original wildtype values.
The TRAF6 KO simulations matched the experimental data, i.e., the removal of TRAF6 only slightly downregulated NF-kB, JNK and p38 activation ( Figure 3B-D TRAF6 KO). For TRADD KO, NF-kB activation was slightly increased in experiments ( Figure 1A). This result was recapitulated by our model which suggests the enhancement of NF-kB activation in TRADD KO is due to signaling flux redistribution [18] (Figure 3B TRADD KO). However, in contrasts to experimental results, which showed only a small downregulation, simulations of JNK and p38 showed almost abolished activation ( Figure 3C-D TRADD KO). Thus, our simulation failure suggests, in actual cells, a novel pathway might exist that compensate the loss of MAP Kinases activation in TRADD KO.
Macrophages with TRIF mutation treated with poly (I:C) showed abolishment of MAP kinase ERK activation [32] and all MAP kinases showed similar kinetics in wildtype and TRAF6 KO [11]. Thus, we added a pathway from TRIF to MAP kinases into our model ( Figure 4A and Table 1) and adjusted the parameter values between TRIF to RIP1, TRADD and the novel pathway to fit wildtype NF-kB, JNK and p38 activation. We re-performed the simulations of TRAF6 KO and TRADD KO. The final model simulations matched experimental outcome for all three conditions ( Figure 4B-D).
To investigate other possible pathways for JNK and p38 activation in poly (I:C) stimulation, reports indicate two other receptors known to recognize dsRNA; the retinoic-acid-inducible protein (RIG)-I and melanoma-differentiation-associated gene (MDA) 5 and, therefore, might be potential candidates  [33,34,35]. However, firstly, RIG-I is unable to recognize poly (I:C) [33,35] and secondly, MDA5 signaling pathway is unknown to trigger MAP kinases [35]. Furthermore, TRIF mutation showed abolishment of ERK activation [32] and all MAP kinases seem to possess similar kinetics in poly (I:C) stimulation [11]. Taken together, these data suggest the predicted pathway for JNK and p38 activation is through TRIF and not by RIG-1 or MDA5. To find the possible candidates to be involved in the novel pathway, we again surveyed the literature. The TRAF family members, six to-date, are well-known to bind to the TIR domain of TRIF with their C-terminal [36,37,38,39], while two members of the RIP family, RIP1 and RIP3, are found to interact with TRIF through the RHIM domain found in RIP1, RIP3 and TRIF [40,41]. RIP1, TRAF6 and TRAF3 already exist in the current TLR3 signaling, while TRAF1 was recently found to inhibit the TRIF-mediated signaling [38,39]. Thus, we suggest that RIP3, TRAF2, TRAF4 or TRAF5 may be part of the novel pathway activating JNK and p38.
In summary, the analysis of the poly (I:C) stimulation in macrophages reveals the involvement of uncharacterized missing intermediary steps prior to TLR3 activation and the existence of a key pathway from TRIF to JNK and p38, but not NF-kB activation. Although further experimental validation is required to identify the uncharacterized processes and molecules nevertheless, we show that the simple mass-action linear response can represent activation dynamics of the TLR3 signaling, and can be used to predict novel features of complicated immune pathways.

Development of in silico TLR3 model using perturbationresponse approach
The basic principle behind our perturbation-response approach is to induce a controlled perturbation of a certain (input) reaction molecule of a system, which is kept at steady-state, and monitor the response of the concentration/activation levels of other molecules (output) of the system. By finding a relationship between the input perturbation (e.g. poly (I:C) stimulation) and output response (e.g. JNK, NF-kB activation), the mechanistic properties of the system such as network topology can be uncovered (e.g. novel TRIF-MAP Kinases pathway) without knowing each individual's detailed kinetics.
To illustrate our approach, let us perturb a stable biological network consisting of n molecules from reference (stable) steadystate. The deterministic kinetic evolution equation is dX i dt~F i X 1 ,X 2 ,::,X n ð Þ ,i~1,::,n, ð1Þ where the corresponding vector form of Eq. 1 is dX dt~F X ð Þ. F is a vector of any non-linear function of the molecules vectors X = (X 1 , X 2 , .., X n ), which represents activated concentration levels of signaling molecules (for example through phosphorylation, binding concentration of transcription factors to promoter regions etc.). The response to perturbation can be written by X = X 0 +dX, where X 0 is reference steady-state vector and dX is relative response from the steady-states (dX t = 0 = 0).
The use of small perturbation around steady-state leads to important simplification to the evolution equation, which can be highly non-linear (Eq. 1), resulting in the approximation of the first-order term. In vector form, ddX dt % LF X ð Þ LX X~X0 dX, noting that zeroth order term F(X 0 ) = 0 at the steady-state X 0 and the Jacobian matrix or linear stability matrix, J~L F X ð Þ LX

X~X0
. The elements of J are chosen by fitting dX with corresponding experimental profiles and knowing the network topology (e.g., activation causality). Hence, the amount of response (flux propagated) along a signaling pathway can be determined using the law of mass flow conservation with first order mass-action kinetics, ddX dt~J dX (Table 1) [13,18,42]. We considered all signaling response to be linear and sequential up to 90 min after poly (I:C) stimulation (equation [2]) as apparent from the formationdepletion activation profiles (see Experiments, and Results and Discussion). Although there will not be any issue for parameter sensitivity in our system because of linear events, multiple solutions of parameter space, Js, can occur [13]. To overcome this, we compare our (wildtype) model with two other (TRAF6 KO and TRADD KO) experimental conditions (see Verification of Model Parameters).

Modeling strategy
We first curated the murine macrophage TLR3 signaling topology from current literature/database sources and by checking with relevant experimental data ( Figure 5, step 1 and see Results and Discussion) [11,12,31,24,40,41,43]. Next, using the topology, we developed a dynamic model of TLR3 signaling on the E-Cell simulation platform [44] to simulate the dynamics of NF-kB and JNK activation ( Figure 5, step 2). The parameters of our model (elements of Jacobian matrix J) were estimated by fitting the simulation profiles of NF-kB, JNK and p38 with corresponding activation profiles semi-quantified from immunoblots, obtained from wildtype macrophages [11] (Figure 5, step 3) (Note that our model works quantitatively if appropriate data is available). If the model successfully fits the wildtype profile of NF-kB and JNK, we accept the model and call it the Reference Model ( Figure 5, step  4). Otherwise, if the model fails, we adjust the parameters values or model topology until the model is able to simultaneously fit the experimental profiles of NF-kB, JNK and p38 reasonably ( Figure 5, step 5).
The model simulation results were deemed acceptable by comparing the simulation profile with the experimentally semiquantitative profile based on three criteria; i) time of activation onset, ii) time to reach peak, and iii) decay profile (Figure 5 insert). For example, JNK experimental profile shows null activation until 10 min, peak activation at 30 min and decaying trend until 90 min ( Figure 1B, WT). In our first simulation of JNK ( Figure 2C, WT), the time of activation onset was almost instantaneous, peak at around 10 min, and complete decay after 60 min. This simulation, therefore, is considered BAD based on the three criteria which poorly matched ( Figure 5 insert, upper panel). On the other hand, after modifying the signaling topology by adding novel intermediary steps, the simulation improved on all three criteria and considered GOOD ( Figure 3C, WT, and Figure 5 insert, lower panel).

Verification of Model Parameters
There can be several parameter spaces for Jacobian matrix, Js, which could fit wildtype semi-quantitative activation profiles of NF-kB, JNK and p38. To reduce the parameter spaces, we performed an iterative two-mode process of model ''creation'' and ''testing''; ''Creation'' is to develop and fit our model with published available semi-quantitative profiles of NF-kB and MAP kinases activity in wildtype (step 4). ''Testing'' is to optimize reaction networks and parameter values, obtained from step 4, by simulating the knockout (KO) conditions (TRAF6 KO and TRADD KO) and comparing the results with their respective experimental profiles of NF-kB and MAP kinases activation [11,12]. Each KO model was generated from the wildtype model by setting the parameters involving the KO molecule null. If our KO simulation does not compare well with experiments, we retune the model parameters till the model simulations fit the activation profiles of both wildtype and TRAF6 KO ( Figure 5 step  6). If unsuccessful, we modify the topology so that a better fit can be obtained, for example, adding novel signaling intermediates to obtain the delayed activation [17] or crosstalk mechanism to provide an alternative source of activation in the KO condition [42,45,46] (Figure 5 step 4, see Results and Discussion). To avoid false prediction, we further check comparing the model simulation against TRADD KO. When the model fails at these steps, it implies that either J (parameter space) is not correct or the represented topology is insufficient. We repeat the procedure of modifying parameter values and topology until we obtained reasonable predictions for all 3 conditions (wildtype, TRAF6 KO and TRADD KO) at 6 time points (0, 10, 20, 30, 60 and 90 min). At this point, we accepted the model and called it the Robust Model ( Figure 5 step 7). Therefore, using steps 1-7, we were able to identify missing or incorrect feature(s) of the TLR3 pathways. The final (robust) TLR3 model's equations and parameter values can be found in Table 1.

Experiments on macrophages
We utilized our published time-course experimental data of NF-kB, JNK and p38 of wildtype (model creation) and TRAF6deficient (model testing) macrophages with 10 mg/ml poly (I:C) stimulation [11]. For further model testing, we utilized published experimental profile of NF-kB, JNK and p38 of TRADD-deficient macrophages with similar poly (I:C) stimulation [12]. The activation levels of NF-kB, JNK and p38 ( Figure 1A-B, Figure  S1A) were quantified from the western blots using ImageJ (http:// rsbweb.nih.gov/ij/).