Mathematical Modeling of Pro-and Anti-Inflammatory Signaling in Macrophages

Inflammation is a beneficial mechanism that is usually triggered by injury or infection and is designed to return the body to homeostasis. However, uncontrolled or sustained inflammation can be deleterious and has been shown to be involved in the etiology of several diseases, including inflammatory bowel disorder and asthma. Therefore, effective anti-inflammatory signaling is important in the maintenance of homeostasis in the body. However, the inter-play between proand anti-inflammatory signaling is not fully understood. In the present study, we develop a mathematical model to describe integrated proand anti-inflammatory signaling in macrophages. The model incorporates the feedback effects of de novo synthesized pro-inflammatory (tumor necrosis factor α; TNF-α) and anti-inflammatory (interleukin-10; IL-10) cytokines on the activation of the transcription factor nuclear factor κB (NF-κB) under continuous lipopolysaccharide (LPS) stimulation (mimicking bacterial infection). In the model, IL-10 upregulates its own production (positive feedback) and also downregulates TNF-α production through NF-κB (negative feedback). In addition, TNF-α upregulates its own production through NF-κB (positive feedback). Eight model parameters are selected for estimation involving sensitivity analysis and clustering techniques. We validate the mathematical model predictions by measuring phosphorylated OPEN ACCESS


Introduction
Inflammation is a beneficial self-defense mechanism that is initiated by the body to eliminate pathogens and prevent the spread of infection [1].The inflammatory responses to pathogens and other inflammatory stimuli are mediated by innate (dendritic cells and macrophages) and adaptive immune cells (T-cells and B-cells) [2].Immune cells have transmembrane receptors called Toll-like receptors (TLR) that recognize foreign molecules based on pathogen-associated molecular patterns (PAMPs), such as flagellin of bacterial flagella [3], lipopolysaccharide (LPS) of Gram-negative bacteria and peptidoglycan of Gram-positive bacteria [4].Recognition of PAMPs by immune cells (such as macrophages) triggers the production and secretion of pro-inflammatory cytokines, which leads to the recruitment of phagocytic cells, such as neutrophils [5], for eliminating pathogens.While inflammation is a beneficial body response, unabated (chronic) inflammation is deleterious, as it can result in immune cells attacking other host cells.Chronic inflammation has been shown to be involved in the etiology of several diseases, including inflammatory bowel disease (IBD) [6] and asthma [7].Chronic inflammation can also arise in the absence of pathogen infection.Since the mucosal immune cells in the gastro-intestinal (GI) tract are in close proximity with intestinal microbiota [8], any alteration in the intestinal microbial community (i.e., dysbiosis) can also lead to uncontrolled pro-inflammatory responses.This sustained inflammation in the absence of any infection has been shown to result in ulcerative colitis or Crohn's disease [6].
Nuclear factor-κB (NF-κB) is an important transcription factor that plays a pivotal role in mediating inflammatory responses in immune cells, such as macrophages [9].NF-κB is made up of two subunits, p50 and p65 [10], and is sequestered as an inactive complex in the cytosol by an inhibitor protein, IκBα [9].When macrophages detect the presence of bacteria (by detecting LPS) through their cell surface receptor, TLR4, an LPS-TLR4 complex is formed that triggers the activation of IκBα kinase (IKK), resulting in phosphorylation of IκBα-NFκB and subsequent ubiquitination and degradation of IκBα [9].NF-κB, which is catalytically released from the inactive IκBα-NFκB complex, translocates into the nucleus and binds to response elements in the promoter region of its target genes to activate their transcription [9].Several target genes with functions in inflammation and immune regulation have been identified for NF-κB [11], of which TNF-α and IL-10 are the most prominent pro-and anti-inflammatory cytokines, respectively [12][13][14].In addition to TNF-α and IL-10, other NF-κB responsive genes that have significant NF-κB regulatory functions are IκBα (sequesters free NF-κB) [15] and A20 (inactivates IKK) [16].However, NF-κB is not the only transcription factor that regulates IL-10 and TNF-α signaling and often acts in concert with other transcription factors.For example, signal transducer and activator of transcription 3 (STAT3) is a well-studied transcription factor involved in IL-10 signaling [17,18].
STAT3 not only regulates transcription of IL-10, but is itself activated by IL-10 [19] and LPS [20] in a feedback manner.The effects of pro-inflammatory cytokines, such as TNF-α and IL-1β, are countered by signaling initiated by anti-inflammatory cytokines.IL-10 is a potent anti-inflammatory cytokine and suppresses the production of pro-inflammatory cytokines, like TNF-α [21], by downregulating NF-κB through inhibition of IKK activation and suppression of free phosphorylated NF-κB translocation from cytosol to nucleus [22,23].Several computational models of inflammatory signaling have been previously developed.These include a model for the IL-6 signal transduction pathway by Singh et al. [24], the TNF-α signaling pathway by Huang et al. [25], Lipniacki et al. [26], Rangamani et al. [27] and Hoffmann et al. [28].A characteristic feature of these models is that they describe the dynamics of signaling initiated by a single pro-inflammatory cytokine.Moya et al. [29] developed a mathematical model to represent interactions between IL-6 (pro-inflammatory) and IL-10 (anti-inflammatory) in hepatocytes when both of these cytokines were used as stimuli to the cells.The current work describes an interplay between de novo synthesized pro-inflammatory (TNF-α) and anti-inflammatory (IL-10) cytokines in macrophages exposed to LPS (Figure 1).Since the inter-play between the pro-and anti-inflammatory signaling in macrophages is poorly understood, our integrated model represents a first step towards modeling the interaction between pro-and anti-inflammatory signaling mediators that is important in inflammation and maintaining homeostasis.

Model Formulation
The mathematical model presented in this paper is an integration of an inflammatory module and an anti-inflammatory module.The model is developed by representing biochemical reactions involved in the signal transduction pathway (Figure 2) as a set of non-linear ordinary differential equations (ODE) of the form: where x is a vector of states, u is a vector of inputs and p is a vector of parameters.The model comprises 29 differential equations (Table 1) and 37 parameters (Table 2).Each differential equation represents the rate of change of the concentration of a particular protein involved in the pathway.[ The inflammatory (TNF-α) module is adapted from Huang et al. [25] and Lipniacki et al. [26].While these models use TNF-α as the input, our model describes LPS (input)-induced signaling through TLR4 (LPS receptor), which leads to TNF-α production.Besides adding TLR4 to the model, the TNF-α receptor description is retained to represent the positive feedback of de novo synthesized TNF-α on NF-κB regulation.We have included a kinetic term for TNF-α mRNA transcription initiated by nuclear NF-κB and component balances for TNF-α in the cytoplasm and the supernatant.In addition, Lipniacki et al., included TRADD, TRAF2, RIP-1, FADD, caspase-3 and caspase-8 proteins, which are left out of the model presented here, as we focused only on some of the key biochemical reactions involved in LPS-induced NF-κB activation, its effect on the production of TNF-α and IL-10 and, in turn, the role of these cytokines on the feedback regulation of NF-κB.The similarities between the model described in Lipniacki et al. [26], and our current ODE model lie in the formulation of the biochemical reactions involved in IKK activation, IκBα-NFκB phosphorylation, dissociation and nuclear transport of NF-κB, nuclear NF-κB-induced IκBα, A20 mRNA transcription, free NF-κB sequestration by de novo synthesized IκBα and IKK inactivation by A20.We added a balance for phosphorylated IκBα, as it is known to degrade after dissociation from the IκBα-NFκB complex.
The anti-inflammatory (IL-10) module is adapted from the IL-6 and IL-10 model by Moya et al. [29].Only the ODEs involved in IL-10 signaling through the IL-10 receptor (as mentioned in Moya et al. [29]) are included in the anti-inflammatory module of our current model to formulate the feedback effects of IL-10 on its own production (through positive feedback regulation of STAT3) and TNF-α production (through negative feedback regulation of NF-κB).Biochemical reactions, as described in Moya et al., for STAT3 phosphorylation, dimerization and nuclear translocation to initiate transcription are retained in our current model.Transcription and translation of SOCS3 due to STAT3 and downstream biochemical reactions associated with SOCS3 are not included in the model presented here.A Michaelis-Menten-type kinetics for IL-10 transcription, initiated by the transcription factors, NF-κB and STAT3, and component balances for IL-10 in the cytoplasm and supernatant, have been included here.
Some values of the parameters (Table 2) and initial concentrations of proteins (Table 3) are adapted from the TNF-α signaling models by Huang et al. [25], Lipniacki et al. [26], Rangamani et al. [27], Hoffmann et al. [28] and the IL-6 and IL-10 model by Moya et al. [29].The previously developed models consisted of 37 differential equations and 60 parameters for the TNF-α model by Huang et al. and 68 differential equations and 118 parameters for the IL-6 and IL-10 model by Moya et al. [29] Among the proteins included in these models, very few are quantifiable by experimental methods, making parameter estimation difficult.In our current integrated model, we have reduced the number of differential equations to 29 and the number of parameters to 37 by only focusing on the key proteins of the pathway.Using a smaller model increased parameter identifiability and simplified parameter estimation.
The ODE model is structurally divided into pro-inflammatory (TNF-α) and anti-inflammatory (IL-10) modules that are both initiated by LPS stimulation and NF-κB activation.Below is the description of the implemented reaction network as shown in Figure 2. Pro-Inflammatory Module: (1) Exogenous LPS binds to the cell surface receptor (TLR4); (2) LPS-TLR4 complex initiates activation of IKK neutral to IKK active ; (3) IKKactive phosphorylates IκBα-NFκB and initiates dissociation of the inactive IκBα-NFκB complex into phosphorylated IκBα and NF-κB species; (4) Free phosphorylated IκBα undergoes ubiquitination and degradation, whereas free cytoplasmic NF-κB translocates into nucleus; (5) Nuclear NF-κB binds to response elements in the promoter regions of TNF-α, IκBα and A20 genes and leads to transcription and translation of the corresponding proteins and subsequent secretion of TNF-α into the supernatant.de novo intracellular IκBα sequesters both free cytoplasmic and nuclear NF-κB by binding to them, and A20 catalyzes the change of IKK active to the IKK inactive form; (6) Secreted TNF-α in the cell culture supernatant binds to its cell surface receptor to form a complex that initiates similar pathways as LPS, resulting in the production of more TNF-α through a positive feedback regulation on NF-κB.
Anti-Inflammatory Module: (1) Nuclear NF-κB binds to the IL-10 gene promoter and initiates transcription of IL-10 mRNA and subsequent translation into IL-10 protein in the cytoplasm, which gets secreted into the supernatant; (2) IL-10 secreted into the supernatant binds to its cell surface receptor, forming a ligand-receptor complex that inhibits activation of IKK neutral to IKK active and translocation of activated free cytoplasmic NF-κB into the nucleus, as well; (3) The IL-10 + receptor complex activates a second transcription factor, presumably STAT3, which, in turn, regulates transcription of the IL-10 gene in a feed-forward manner.Different LPS concentrations (0, 0.1, 1, 10 μg/mL) are used to stimulate the model.The IKK complex and NF-κB dimer (p50-p65) are considered as single proteins in the model.The signal transduction model comprises feedback regulatory loops involving TNF-α and IL-10.The positive feedback of TNF-α on its own production is represented by de novo TNF-α binding to its cell surface receptor, activating IKK neutral to IKK active , leading to phosphorylation and dissociation of the IκBα-NFκB complex to release NF-κB, which translocates to the nucleus to initiate transcription of TNF-α.IL-10 has a negative feedback effect on TNF-α production by inhibiting NF-κB activation (phosphorylation and dissociation), and the extent of this inhibition is calculated on the basis of the ligand bound IL-10 receptor complex (IL10-IL10R) concentration, which is represented as kin (Equation ( 2)): The maximum attainable concentration of IL10-IL10R is denoted by IL10-IL10Rmax with an assumed value of 2.56 × 10 −6 .kin is multiplied by factors that are inhibited by IL-10, such as IKK activation and nuclear translocation of free cytoplasmic NF-κB [22,23] (Table 1).The higher the concentration of the IL10-IL10R complex, the lower will be the value of kin and, hence, the lower will be the contribution of the terms mentioned above to the total outcome of NF-κB signaling, resulting in suppression of TNF-α production by IL-10.Positive feedback of IL-10 on its own production is represented by a set of differential equations that describe the IL-10 bound receptor complex phosphorylating transcription factor STAT3, which then dimerizes and translocates into the nucleus, binds to the promoter region of the IL-10 gene and initiates transcription of IL-10 mRNA and subsequent translation and secretion of IL-10 protein.

Parameter Selection and Estimation
The parameter estimation problem for a dynamic system described by ordinary differential equations (ODEs) can be mathematically formulated as follows: Subject to, ( ) = ( , , ), (0) = (4) ≤ ≤ ≤ ≤ (7) where yik and ŷik are the simulated and measured output data of the i-th component at sampling time tk, respectively (Equation ( 3)); p are the parameters to be estimated, which are selected by local sensitivity analysis; x are the state variables of the dynamic system with initial values x0; and u are the inputs to the system (Equation ( 4)).In addition, the state variables x and parameters p are restricted within certain ranges, as shown in Equations ( 6) and ( 7), determined by the underlying biology and prior knowledge based on mathematical models developed by Lipniacki et al. [26], Huang et al. [25] and Rangamani et al. [27].The simulated output vector y (Equation ( 5)), which is validated by experimental data, includes the intracellular ratio of phosphorylated NF-κB to total NF-κB (relative to the control) and the concentration of the cytokines, TNF-α and IL-10, in the cell culture supernatant.Since the experiments are conducted with four different levels of the input u (i.e., different concentrations of LPS), four sets of measured outputs ŷik are obtained.Three sets of data obtained for 0, 0.1 and 1 µg/mL LPS stimulations are used for parameter estimation, and the fourth dataset for 10 µg/mL LPS stimulation is used for model validation.
First, the set of parameters that are to be estimated are selected by local sensitivity analysis and hierarchical clustering.Following this, the trust-region optimization technique is used to estimate the selected parameters.This technique is used in this work, as it is able to handle singular Hessian matrices, significant uncertainty in the parameters of the models and noisy data.In this technique, an optimization algorithm is applied in an outer loop, while the evaluation of the objective function and its gradients are performed by numerical integration of the ODEs in the inner loop [32,33] (shown in Figure 3).The trust-region method is guaranteed to converge to local optima with much weaker assumptions than line search methods.In this work, fmincon (MATLAB function) is used as the NLP (non-linear programming) solver and ode15s (MATLAB function) is used as the ODE integrator.It is worth noting that ode15s is specifically designed for stiff systems, such as the model discussed here, where both fast and slow dynamics exist, e.g., in our system, phosphorylation of NF-κB is a much faster process than de novo synthesis of TNF-α and IL-10.

Cell Culture and Experimental Set-Up
The murine macrophage cell line RAW264.7 (gift from Dr. Paul deFigueiredo, Texas A & M University) was routinely cultured in DMEM with 10% FBS.LPS (heat-killed Salmonella enterica) was purchased from Sigma Aldrich (St. Louis, MO, USA).

Transcription Factor NF-κB Quantification by ELISA
The relative concentrations of total and phosphorylated NF-κB in LPS-stimulated RAW264.7 macrophage cells were determined using a commercially-available enzyme-linked-immunosorbent assay (ELISA) kit (R&D Systems, Minneapolis, MN, USA) according to the manufacturer's suggested protocol.

Cytokines TNF-α and IL-10 Quantification by ELISA
The concentrations of de novo synthesized TNF-α and IL-10 in LPS-stimulated RAW264.7 macrophage culture supernatant were determined by commercially available enzyme-linked immunosorbent assay (ELISA) kits (Thermo Scientific, Rockford, IL, USA), using the manufacturer's suggested protocol.

Results and Discussion
Based on published reports of LPS stimulation resulting in TNF-α [34] and IL-10 secretion [35] in RAW264.7 murine macrophages, as well as the established suppression of TNF-α by IL-10 in RAW264.7 cells [34], we developed an integrated ODE model to represent the production of TNF-α and IL-10 in RAW264.7 cells upon LPS stimulation and their regulatory feedback loops.Eight parameters of the model are selected for estimation using local sensitivity analysis and hierarchical clustering (shown in Figure 4) [36,37].The y-axis represents the parameter distance ranging from zero to one (the larger the distance, the smaller the similarity between the parameters).The red line presents the cutoff value, which groups the entire set of parameters into eight pairwise indistinguishable clusters.The selected parameters, which have the largest sensitivity magnitude in each cluster, are highlighted in red.The values of the selected parameters are estimated using the trust-region optimization technique, as described in Materials and Methods section, and their estimated values are listed (in bold) in Table 2.One advantage of this approach is that the selected parameters used for estimation result in a more robust dynamic model with an accurate prediction capability [37].Model simulations after parameter estimation predict rapid phosphorylation of NF-κB upon exposure to LPS, as shown in Figure 5A.The maximum fold change in the ratio of phosphorylated NF-κB to total NF-κB in LPS-treated cells relative to control increased with increasing LPS concentration and varied from ~2.0 at the highest LPS concentration to being essentially unchanged at the lowest concentration (Figure 5A).Simulation of de novo synthesized TNF-α profile upon LPS stimulation shows that the TNF-α concentration reaches a maximum of ~1500 pg/mL at 4 h and starts declining thereafter.Even though LPS is continuously present, TNF-α is undetectable at 24 h for all LPS concentrations (Figure 5B).The maximum TNF-α concentration at 4 h increases with increasing concentrations of LPS.According to the model predictions, the de novo synthesized IL-10 concentration increases beyond 2 h of LPS stimulation, as shown in Figure 5C, and the concentration of IL-10 produced increases with increasing LPS concentrations.Model simulations are validated by experimental data obtained from LPS-stimulated RAW264.7 cells.The comparison between simulated and experimental data for the fold change in NF-κB (phosphorylated NF-κB/Total NF-κB, relative to control), de novo TNF-α and IL-10 concentration profiles after parameter estimation are shown in Figure 5A, 5B and 5C, respectively.In Figure 5A, the discrepancy between the model simulation and experimental data for 0.1 µg/mL LPS stimulation could arise from our assumption that the binding affinity between LPS and TLR4 is constant and concentrations of the LPS-TLR4 complex are linearly proportional to the concentrations of LPS tested.However, in reality, the ligand-receptor binding kinetics might follow a non-linear behavior, which is not accommodated in our computational model.The binding of LPS to TLR4 even at lower concentrations of LPS (e.g., 0.1 µg/mL) might result in higher concentrations of the LPS-TLR4 complex, resulting in more downstream phosphorylation of NF-κB (as indicated by the higher phosphorylated NF-κB/total NF-κB ratio for the experimental data in Figure 5A) than the model is able to predict.However, the dynamic model where the parameters have been estimated exhibits a reasonably good fit for phosphorylated NF-κB/total NF-κB profiles (relative to control) at 0 μg/mL, 1 μg/mL and 10 μg/mL LPS stimulations.Furthermore, the model exhibits reasonable agreement between simulated and experimental data for both training and validation datasets for the TNF-α and IL-10 dynamic concentration profiles.
The model suggests that the initial increases in TNF-α and IL-10 are due to NF-κB activation (i.e., phosphorylation and dissociation of NF-κB from IκBα-NF-κB complex) in the cytoplasm and subsequent gene expression in the nucleus; and the decrease in TNF-α concentration after 4 h is due to negative feedback of IL-10 on NF-κB activity (by inhibiting both IKK activation and nuclear translocation of phosphorylated NF-κB).Interestingly, the levels of IL-10 continue to increase, even when the levels of activated NF-κB are no longer increasing.It is possible that the initial burst of IL-10 produced through NF-κB activation can activate other transcription factors, which leads to an increase in IL-10 levels.An example of this could be the transcription factor STAT3, which has been shown to be activated by IL-10 [38].Endogenous IL-10 in LPS stimulation of macrophages (RAW264.7)[39] forms the IL-10-IL10R complex that initiates phosphorylation of cytosolic STAT3, followed by its dimerization and translocation into the nucleus.The STAT3 dimer binds to the DNA response element and triggers transcription of IL-10 (Figure S1).Thus, increasing LPS concentration could lead to increasing concentrations of IL-10 due to the positive feedback of IL-10 on its own production.
It can be seen that the model predictions and experimental data are in reasonable agreement, which demonstrates that biochemical reactions, which form the structure of the model, are physiologically relevant and can depict the interplay between pro-inflammatory and anti-inflammatory immune responses to maintain equilibrium (homeostasis).Furthermore, cross-talk between positive and negative feedback regulatory loops, incorporated in the model, is integral towards mathematically representing biochemical and gene regulatory networks, as mentioned by Tian et al. [40].Our model also shows that the anti-inflammatory functions in the RAW264.7 macrophage cell line is initially triggered by pro-inflammatory stimulation.This model structure can be extended to study other cell types with a modification in parameter values to fit model predictions to experimental datasets for specific cell types.
The integrated mathematical model of pro-and anti-inflammatory host immune response discussed in this paper is a step towards assimilating our knowledge and developing a quantitative understanding of the signal transduction pathways involved in maintaining immune homeostasis, disruption of which can lead to inflammatory disorders.This mathematical model can be further used to study intra-and inter-kingdom signaling, i.e., the effect of bacterial metabolites (e.g., indole) synthesized by micro flora present in the host gastro-intestinal tract on host immune response [41].This model can be used as the basic structure to incorporate additional transcription factors, which will be needed to study the signaling of indole and their interactions with NF-κB.

Figure 1 .
Figure 1.Schematic representation of NF-κB signal transduction pathway under LPS stimulation in macrophages.LPS binds to cell-surface TLR4, forms the LPS-TLR4 complex that initiates activation of IKK, subsequent rapid phosphorylation and dissociation of the IκBα-NFκB complex.Phosphorylated IκBα undergoes degradation, whereas free cytoplasmic NF-κB translocates into the nucleus, binds to DNA response elements and initiates the transcription of target genes TNF-α, IL-10, IκBα and A20.De novo synthesized TNF-α and IL-10 are secreted into the cell culture supernatant, where they bind to their respective cell surface receptors and initiate their positive (TNF-α) and negative (IL-10) feedback regulations on NF-κB.The LPS-induced NF-κB signaling pathway is indicated in solid blue arrows.TNF-α-induced positive feedback regulation of NF-κB is indicated in dashed cyan arrows, and IL-10-induced negative feedback regulation of NF-κB is indicated in solid red lines.

Figure 3 .
Figure 3. Algorithm for parameter estimation used to optimize model parameters.An optimization algorithm is applied in an outer loop, while the evaluation of the objective function and its gradients are performed by numerical integration of the ODEs in the inner loop.fmincon (MATLAB function) is used as the NLP (non-linear programming) solver and ode15s (MATLAB function) is used as the ODE integrator.

Figure 4 .
Figure 4. Representation of local sensitivity analysis results, used for selecting parameters that are to be estimated.The y-axis represents parameter distance ranging from zero to one.The red line represents the cutoff value, which groups the entire set of parameters into eight pairwise indistinguishable clusters.The selected parameters from each of the eight clusters are highlighted in red.The normalized sensitivity magnitudes of the parameters are reflected in the histograms.

Table 1 .
Differential equations representing biochemical reactions involved in LPS-induced NF-κB signal transduction pathway, as used in the ODE model. 1.

Table 2 .
List of parameters used in the ODE model.

Table 3 .
State variables and their initial values as used in the ODE model.