Skip to main content

Advertisement

We’d like to understand how you use our websites in order to improve them. Register your interest.

Deriving a Boolean dynamics to reveal macrophage activation with in vitro temporal cytokine expression profiles

Abstract

Background

Macrophages show versatile functions in innate immunity, infectious diseases, and progression of cancers and cardiovascular diseases. These versatile functions of macrophages are conducted by different macrophage phenotypes classified as classically activated macrophages and alternatively activated macrophages due to different stimuli in the complex in vivo cytokine environment. Dissecting the regulation of macrophage activations will have a significant impact on disease progression and therapeutic strategy. Mathematical modeling of macrophage activation can improve the understanding of this biological process through quantitative analysis and provide guidance to facilitate future experimental design. However, few results have been reported for a complete model of macrophage activation patterns.

Results

We globally searched and reviewed literature for macrophage activation from PubMed databases and screened the published experimental results. Temporal in vitro macrophage cytokine expression profiles from published results were selected to establish Boolean network models for macrophage activation patterns in response to three different stimuli. A combination of modeling methods including clustering, binarization, linear programming (LP), Boolean function determination, and semi-tensor product was applied to establish Boolean networks to quantify three macrophage activation patterns. The structure of the networks was confirmed based on protein-protein-interaction databases, pathway databases, and published experimental results. Computational predictions of the network evolution were compared against real experimental results to validate the effectiveness of the Boolean network models.

Conclusion

Three macrophage activation core evolution maps were established based on the Boolean networks using Matlab. Cytokine signatures of macrophage activation patterns were identified, providing a possible determination of macrophage activations using extracellular cytokine measurements.

Background

Over the past 30 years, extensive research has been dedicated to investigating the role of macrophages due to its versatile functions in innate immunity, infectious diseases, and progression of cancers and cardiovascular diseases, the top 2 leading causes of death in the world [1,2,3,4]. Macrophages have bactericidal and phagocytic functions and regulate immune responses and the development of inflammation by secreting monokines including enzymes, complement proteins, and regulatory factors such as IL-1 (Interleukin-1), IL-6, and IL-12 [5]. M2 macrophages turn off the damaging immune system activation, function in constructive processes like wound healing and tissue repair, and coordinate the chronic inflammatory response by regulating downstream cellular functions (inflammation resolution, endothelial cell, and fibroblast activation) [4, 6,7,8]. These contradictory functions of macrophages were conducted by different macrophage phenotypes classified as M1 (classically activated macrophages) and M2 (alternatively activated macrophages) due to different stimuli in the complex in vivo cytokine environment.

Macrophages have shown in vitro M1 activation with Lipopolysaccharide (LPS), interferon-gamma (IFN-γ), or TNF (Tissue Necrotic Factor) stimuli. M1 macrophages secrete high levels of IL-12 and IL-23 and low levels of IL-10. In contrast, M2 macrophages secrete high levels of IL-10, Transforming Growth Factor–β (TGF-β) and low levels of IL-12 and IL-23. Recently, M2 activation has been further classified into 4 sub-phenotypes, M2a, M2b, M2c, and M2d. M2a activation is stimulated by interleukin-4 (IL-4) and IL-13, while M2b with immune complex (IC) + Toll-like receptor (TLR), IC + IL-1 receptor (IL1R), or IL-1β, M2c with IL-10 or TGF-β stimuli, and M2d with TLR ligands or adenosine receptor ligands [9,10,11,12]. In addition, switches between M1 and M2 phenotypes have been reported [13,14,15].

Analysis of macrophage activations has been reported by others and our group. We have reported differential equation models on macrophage activation and effects of matrix metalloproteinase-9, and − 28 on macrophage activations [16,17,18]. Martin and colleagues have validated a Boolean dynamics of genetic regulatory network for LPS-stimulated macrophage activation using temporal microarray (mRNA) data at 8-time points (15mins, 30mins, 1 h, 2 h, 4 h, 8 h, 16 h, 32 h) [19]. Rex’s group characterized the inflammatory gene expression patterns of LPS-stimulated M1 activation and IL-4 and IL-13-stimulated M2 activation at 0.5, 1, 2, 6, 10, and 24 h, and modeled macrophage activations by combining Boolean dynamics and differential equations [20]. Though less than 50 genes were considered in each macrophage activation model, the limited 6 temporal measurements still led to an unavoidable difficulty in parameter estimation for differential equation models. Recently, transcriptional and post-transcriptional graphical networks of macrophage in healthy and diseased hearts were reported [21]. All these models focused on the regulatory networks at the gene level and gave a deeper insight of what happened inside a cell. Interestingly, all models reported significant cytokine markers for macrophage activations. Furthermore, in an in vivo situation, stimuli and final secretion products of activated macrophages stay outside the cells. It’s natural and practical to define macrophage phenotypes with cytokine expressions for an Input-Output representation only based on cytokine expressions outside of the cells. Therefore, in vitro cytokine profiles from macrophage activations serve perfectly to explore such abstract models.

We globally searched “temporal” or “time series” and “macrophage activation” or “macrophage polarization” in PubMed, GEO datasets/profiles, and ProteomeExchange to screen specific macrophage activation temporal profiles with different stimuli. A total of 173 datasets were found for Mus musculus and 156 for Homo sapiens as of June 18, 2019. Eleven studies focused on temporal macrophage activations. Currently, complete temporal in vivo genome or proteome expression profiles of macrophage activation has not been reported. However, in vitro temporal data has been deposited in public databases to shed a light on macrophage activation mechanisms. A dataset from Melton’s group was selected to establish a Boolean model because macrophage activation expression profiles of 27 cytokines for 4 groups (a control group and 3 activation groups with 3 different stimuli, LPS, IL-4, and IL-10) at 7 points (0 h, 0.5 h, 1 h, 3 h, 6 h, 12 h, 24 h) were documented [22]. This dataset provides extra information on time points, sub-activation groups, and cytokine expression levels to existing models on gene regulations [19,20,21].

In this study, expression levels of the selected differentially expressed cytokines were binarized with “1” (up-regulated) and “0” (down-regulated). Linear programming algorithm was applied to determine the structure (links and interaction strengths) of a regulatory network for macrophage activations among the 27 selected cytokines under 3 different stimuli. Based on the obtained network structure, Boolean functions were determined for each cytokine in the network using a Karnaugh Map. For the first time, core networks for macrophage activation models were generated and mathematically represented with the semi-tensor product. The core network for M1, M2a, and M2c only contained 6, 7, and 5 proteins, respectively. Links of macrophage activation networks were validated based on public Protein-Protein Interaction (PPI) and Kegg pathway databases. The core network for each activation pattern was validated by the coincidence of predicted and temporal experimental results. A novel temporal evolution map for each macrophage activation model was also generated to show all Boolean states, suggesting possible evolutionary paths of the networks with different initial conditions.

Results

Curve fitting with smooth spline algorithm

An illustration of 3 typical expression profiles by curve fitting was shown in Fig. 1. The 3 proteins are C-C motif chemokines ligand 12 (CCL12, also named as MCP-5) in M2c activation, CCL7 (MCP-3) in M1 activation, and chemokine (C-X-C motif) ligand 10 (CXCL10, also named as IP-10) in M2a activation. Interestingly, MCP-5 followed a similar profile as promotion, (\( c\frac{t}{1+t}\Big) \), where c represents a constant of its max amplitude with respect to time t (Fig. 1a). MCP-3 followed a bell shape curve similar to a Gaussian distribution in Fig. 1b, and IP-10 followed an inhibition relatable to, (\( c\frac{1}{1+t}+ bias\Big) \), in Fig. 1c. The bias denotes an offset term since an expression may not necessarily reach 0. The average R-Squared error of the fitting algorithm was 0.9087, 0.8566, and 0.8106 for M1, M2a, and M2c group, respectively.

Fig. 1
figure1

Three typical temporal protein expression profiles for macrophage activations. Horizontal coordinate represents time and the vertical coordinate denotes expression levels of proteins. The red line shows the constructed continuous profile while the blue stars denote the real experimental data at each sampling time point. a The profile of MCP-5 in M2c activation shows a promotion pattern. b The profile of MCP-3 in M1 activation shows a bell curve. c The profile of IP-10 in M2a activation shows an inhibition pattern

A temporal profile is considered as a promotion if its final expression level is more than twice or inhibition if it is less than half of the initial expression level, and a bell-curve otherwise. It’s worth to mention that an expression profile demonstrates different patterns with different time periods due to activation dynamics. About 138 out of 324 (27x4x3, 27 protein expression measured from 4 replicates in 3 activation groups) protein expression profiles followed a promotion pattern, 27 had inhibition profiles, and 148 had a bell curve within the first 6 h as shown in Table 1. Examining the complete 24 h time span, 78 profiles fit a promotion, 110 fit an inhibition while 125 for a bell curve, suggesting significant protein up-regulation in the first 6 h and gradually evolve to stable polarizations. M2 activation profiles demonstrate more inhibition patterns compared to M1 activation in Table 1. Chemokine (C-X-C motif) ligand 1 (CXCL1/KC-GRO) did not or barely reached detectable levels in both M2a and M2c activations. These profiles were denoted as N/A in Table 1. Therefore, KC-GRO was not considered in either M2 activation models for undetectable expression levels. GM-CSF (Granulocyte-macrophage colony-stimulating factor) from some replicates reached a detectable level 12 h post-stimulation in M2a and M2c and was included in both M2 models. KC-GRO was significantly up-regulated in M1 activation, illustrating their role to identify M1 activation [23]. All other protein expression levels are either elevated or down-regulated comparing against the average of the control group at least 1 out of 7 time points with p-value < 0.05.

Table 1 Properties of Expression Profiles

Binarization

Both normal K-means and iterative K-means algorithms were applied to the 50 re-sampled data from the fitted curves. Figure 2 gives an illustration of 2 clusters with the normal K-means method and the final binarization of the iterative K-means method. The Iterative K-means generates better binarization results than the normal K-means method since it identifies more significant peaks in the conversion.

Fig. 2
figure2

Comparison of binarization threshold using K-means and iterative K-means methods. Analyzing the same expression profiles in (a), the K-means method identified the first 5 peaks with relatively high expression levels in (b) while the iterative K-means method identified all peaks in (c)

A binarization obtained from the smooth spline and iterative K-means is shown in Fig. 3. The blue stars represent the initial time points (0 h, 0.5 h, 1 h, 3 h, 6 h, 12 h, 24 h) in the dataset, the red solid line represents the curve fitting result, the green solid line gives the threshold from the iterative K-means, and the yellow solid line is the final binarization result. An error band centered at the binarization threshold is denoted by dashed lines for possible binarization errors.

Fig. 3
figure3

Binarization of a protein expression level. The blue stars show the experimental measurements, the red curve is the fitted continuous profile, the green line represents the threshold from the iterative k means method, and the orange line denotes the final binarization. An illustration of a possible sensitivity region is shown between the two dashed black lines as an error band

Statistical analysis

All experimental measurements at 7 sample times have been analyzed in Melton’s paper for statistical significance [22]. We confirmed their analysis with the Kolmogorov-Smirnov test (p < 0.01). We fitted temporal expression profiles for all genes with detectable expression levels in 3 stimulated groups and 1 control group. Two-sample Kolmogorov-Smirnov tests were performed on fitted temporal profiles of each gene given a stimulus versus the control with a significance value p < 0.01.

The fitted smooth spline functions are determined by minimizing the summation of squared errors at seven time points. The average error at these seven sample time points can be obtained in Matlab. Such error may lead to a false determination of binarization value. To test the sensitivity of binarization, an error band was formed as a binarization threshold ± error. If any of the 50 re-sampled data from the averaged fitted expression stays in the error band, there might be a false binarization and the re-sampled data is sensitive to binarization. The maximum likelihood for false binarization is less than 5% for all genes in 3 stimuli groups. No outlier is detected in our experimental and fitted data.

Determining Boolean functions

There are three assumptions in determining the Boolean functions of a network model. 1) Significant functional changes appear in a period larger than the simulation intervals; 2) The cycles of macrophage activation are larger than our simulation interval; 3) The Boolean network will have a fixed structure and regulation during the 24-h response period in the experiment.

Since macrophage activations have the most significant changes in the first 9 h post-stimulation, we put 40 samples during the first 9 h, which led to a simulation interval of about 13.5 min. The other 10 out of 50 samples were evenly re-sampled from 9 to 24 h. Therefore, the simulation interval is 13.5 min in this study.

Currently, most of the published results show a period of 24 h for macrophages’ responses to stimuli. In addition, experimental data are collected with different time periods ranging from 15 min to hours [19,20,21, 24]. Our simulation time interval is very similar to the minimum sampling time of available experimental measurements (15 min). This validated assumptions 1 and 2. Since biological regulations have demonstrated stability and robustness, assumption 3 is a general assumption for modeling.

Boolean functions are determined with respect to binarized values in a time period instead of a single distinct time point. We verified the Boolean function determined by the averaged gene expressions with two methods: agreement with Boolean function dependent on experimental results and agreement of promotion/inhibition determined by Linear Programming. A total of 23 out of 27 datasets (85%) have shown the agreement of Boolean functions determined by the average expressions of genes from 4 replicates with a total of 50 samples for 24 h post-stimuli and linear programming results.

However, variations of temporal responses exist in biological systems [25]. The variations may be more severe during the transition between two different states. If there is any disagreement of Boolean function determined by Karnaugh map and experimental measurement or Linear programming results, we further determined the Boolean function among proteins with respect to each individual experimental replicate instead of the average of 4 replicates and enriched our samples to be 200. An illustrative example showing the agreement of Boolean functions determined by the averaged expression and individual replicates is shown in Fig. 4. In a small network with two inputs (IL-10 and MIP-1β) and one output (MIP-2), the averaged profiles have a missing state (IL-10 = 1 and MIP-1β = 0), so the Boolean function determined by the Karnaugh map showed MIP-2 = 0, contradicting to our experimental results. In this case, expression levels of 4 replicates are considered together to determine the Boolean function. Note that there are only 49 (average) and 196 (4 replicates) sample time points for the iterative evolution ends at sample 49 to predict the state at sample time point 50. The agreement of this Boolean function is about 90%.

Fig. 4
figure4

The Binary Expressions of a simple network with two inputs (MIP-1β and IL10) and one output (MIP-2). The sub-Figs. a-d showed the binarization of 4 replicates of the LPS stimulated responses and e being the binarization of the average of 4 replicates. The green dashed line represents IL-10, the red dotted for MIP-1β and the blue solid line for MIP-2. The Karnaugh maps for the average and each replicate are shown in f and g respectively

Properties of the networks established with linear programming algorithm

Linear programming algorithm was applied to determine the network structure for 3 macrophage activation patterns. As shown in Table 2, the networks obtained were still complex. To simplify the network structure, a link between two nodes was removed if its weight calculated by linear programming was less than 70% of the maximum link weight described in the Method Section. The 70% threshold was chosen based on the criteria to minimize the number of links and maximize the summation of link weights. Properties of the networks established with original linear programming and the simplified network were shown in Table 2. Input nodes are proteins with no parent nodes and output nodes are proteins with no child nodes. All 3 simplified networks have fewer nodes and links than the original results obtained from the linear programming algorithm. In addition, the complexities of M1 and M2a polarization networks are less than the networks for M2c polarization, specifically for the number of links and network density. While the M2c activation network has 3 nodes with ≥4 parent nodes, these 3 nodes are all output nodes and will not affect the regulation of Boolean functions, suggesting that the 50 sample data should be sufficient to establish the Boolean functions.

Table 2 Properties of M1, M2a, and M2c Activation Networks

LPS induced M1 activation network

The structure for the LPS induced M1 activation network was shown in Fig. 5. Boolean regulations for each node were shown in Table 3. An overbar denotes a negation, while a double vertical bar for a logic ‘OR’ and a ‘&’ symbol representing a logic ‘AND’. In addition, interaction strengths calculated by the linear programming algorithm for each link were listed in Table 3.

Fig. 5
figure5

The network for M1 activation includes 26 nodes with each node representing a protein (blue block), a black line with an arrow for a positive connection from a parent node to a child node, and logic negation being a red dashed connection with a ‘T’ arrow. The core network is shown as the diamond nodes in the network

Table 3 Boolean Functions and Description of LPS Stimulated M1 Activation Network. IS: intermediate state

There are 26 proteins in this network. TIMP-1 (Tissue Inhibitor of Metalloproteinases) is the only protein removed from the 27 measured proteins as an orphan node. The two input nodes are FGF-9 (Fibroblast Growth Factor-9) and IL-4. FGF-9 has been known to be an inflammation promoter in multiple sclerosis [26]. Interestingly, IL-4, not LPS, was annotated as input for the M1 activation network. Though IL-4 is known as a stimuli for M2a activation and anti-inflammatory responses, the combination of LPS and IL-4 has been reported for its promotion of other inflammation cytokines including IL-6, CCL-1, CCL-3 (macrophage inflammatory protein 1-α, MIP-1α), CCL-4 (macrophage inflammatory protein 1-β, MIP-1β), CCL-5 (Regulated upon Activation, Normal T cell Expressed, and Secreted, RANTES), TNF-α, and IFN-γ, which are also included in our M1 activation network [44]. In addition, LPS has been shown to induce IL-4 gene expression, which confirmed the rationale for IL-4 serving as an input for the M1 activation network [59]. The M1 activation network has 10 output nodes: GM-CSF, IL-2, IL-6, IL-7, IL-12p70, IL-17A, IP-10, OSM (Oncostatin M), TNF-α, and VEGF-A (Vascular Endothelial Growth Factor A). All these proteins’ expressions are significantly upregulated compared to the control group.

Validation of LPS induced M1 activation network model

To further validate the links in our network model, Protein-Protein-Interaction (PPI) database (STRING), KEGG pathway database, and literature results were used to confirm interactions among proteins in the network. The majority of links (35 out of 39) were confirmed by STRING directly. The other 4 links can be verified through the KEGG pathway database as shown in Table 3. Lymphotactin can trigger the chemokine receptor pathway which can further stimulate multiple pathways such as Salmonella infection to produce IL-1α. MCP-5 to IL-12p70 link can be verified through MCP-5, Chemokine receptor pathway, MAPK signaling pathway, and then RIG-I-like receptor signaling pathway that leads to the synthesis of inflammatory cytokines including Il-12. Similarly, MCP-5 can also trigger the secretion of IL-17 through chemokine receptor pathway, JAK-STAT signaling pathway, and Th17 cell differentiation pathway. MCP-3 to RANTES link can possibly be verified through MCP-5, Chemokine Receptor Pathway, MAPK signaling pathways, Toll-like receptor signaling pathway and Herpes simplex virus 1 infection/ Influenza A pathway.

There were 5 feedback loops as shown in Fig. 5: 1) IFN-γ attractor; 2) from IFN-γ to itself through MCP-1 (monocyte chemoattractant protein), MCP-5, Lymphotactin, IL-1α, and SCF (Stem Cell Factor); 3) from IFN-γ to itself through Lymphotactin, IL-1α, and SCF; 4) a loop between Lymphotactin and MCP-5, and 5) a loop between MCP-3 and RANTES. Interestingly, the first three loops illustrated an IFN-γ attractor self-promotion loop during macrophage activation [60]. STRING database showed MCP-3 and RANTES, MCP-5 and Lymphotactin are interaction pairs.

Figure 5 showed the complete Boolean network generated with our temporal profiles. The 10 output nodes (from the leftmost node IL-7 to the rightmost node VEGF-A at the bottom of Fig. 5) do not affect the regulations of the network and their states can be determined by their parent nodes. In addition, the signal flows from FGF-9 to IL-11 to MIP-1β and from SCF to IL-3 to MIP-1β form single direction forward links, therefore, given the status of input node FGF-9, the status of MIP-1β can be determined together with SCF. Similarly, signal flow from IL-4 to IL-10 to KC-GRO forms single directional links. The input node IL-4 and Lymphotactin expression levels determine IL-10 and KC-GRO expression levels, suggesting Lymphotactin should be a major in evolution and IL-10 and KC-GRO are dependent variables in the simulations and can be removed. Therefore, removing output and single direction transduction nodes led to a simplified network with only 8 proteins for two sub-networks: MCP-3 – RANTES loop for sub-network 1 and IFN-γ, MCP-1, MCP-5, Lymphotactin, IL-1α, and stem cell factor (SCF) for a core 6-protein (in diamond shapes) network as shown in Fig. 5. States of other nodes in the M1 activation network can be determined by the states of these 8 proteins and the input nodes.

The state transition map for the M1 activation core network was shown in Fig. 6. The state transition map contains 64 = 26 states with one isolate state \( {e}_{64}^7 \) and the final state \( {e}_{64}^{58} \). The state \( {e}_{64}^{62} \) is an important transition state with multiple incoming branches. Our binarized and original temporal profiles showed the transition from \( {e}_{64}^{11}\to {e}_{64}^1\to {e}_{64}^{22}\to {e}_{64}^{62}\to {e}_{64}^{58} \) and these states were illustrated in red boxes in Fig. 6. The binary values of these 5 transition states were shown in Table 4. The state numbers are illustrated in two’s complement representation. Temporal transitions for the 6 proteins with binary states and sampled data were shown in Fig. 7.

Fig. 6
figure6

The evolution map of the core network for M1 activation showed all possible pathways predicted by the M1 Boolean network. The evolution of 64 (64 = 26) states for 6 nodes in the core network were illustrated from different initial conditions. The states in red boxes represent the evolution path based on the temporal cytokine profiles in this study

Table 4 Binary Expressions of Proteins in the Core Network for LPS Stimulated Macrophage Activation. The binarization is shown in two’s complement representation
Fig. 7
figure7

Experimental expression levels of the 6 proteins in the core network for M1 activation were compared with their binary profiles. The blue stars showed the measured experimental expression levels at 6 time points (0, 0.5 h, 1 h, 3 h, 6 h, and 12 h) while the orange line illustrated the temporal binary expression. The vertical black bars represent the states that are in our M1 activation evolution map

It’s worth to mention, not all 4 experiments for M1 activation converged at the same time. One replication of this experiment reached the state \( {e}_{64}^{58} \) while the other 3 ended at the state \( {e}_{64}^{62} \). The coincidence of state transition from experimental results and computational results demonstrated the effectiveness of the core network model of M1 activation.

M2 activation networks

Boolean networks for IL-4 induced M2a and IL-10 induced M2c activations have been established and shown in Fig. 8 and Fig. 9, respectively. M2 activations were more complicated than M1 activation due to variant stimuli and M2 subtypes [11, 14]. These were represented by the increased network density and number of links in the M2c activation network compared against the M1 activation network as shown in Table 2.

Fig. 8
figure8

The network for M2a activation includes 24 nodes with each node representing a protein (blue block), a black line with an arrow for a positive connection from a parent node to a child node, and logic negation being a red dashed connection with a ‘T’ arrow. The core network is shown as the diamond nodes in the network

Fig. 9
figure9

The network for M2c activation includes 24 nodes with each node representing a protein (blue block), a black line with an arrow for a positive connection from a parent node to a child node, and logic negation being a red dashed connection with a ‘T’ arrow. The core network is shown as the diamond nodes in the network

Though M2a and M2c networks shared common responses, different subtypes also illustrated specific responses [3]. In the M2a and M2c activation networks, there were 23 common nodes as shown in Tables 5 and 6. The only two different proteins between M2a and M2c are TIMP-1 and IL-11. TIMP-1 was included in both M2a and M2c networks from our linear programming results. However, TIMP-1 was an orphan node in the M2c network and not shown in Table 6 and Fig. 9. Previous results have shown that TIMP-1 is consistently produced by macrophages and further induced by LPS (M1) and IL-10 (M2c) in fully differentiated macrophages [128, 129]. Additionally, TIMP-1 was down-regulated by IL-4 in M2a [130]. Since TIMP-1 functions by forming one-to-one complexes with target metalloproteinases and there is no measurement on metalloproteinases in this study, TIMP-1 is classified as an orphan node. IL-11 is included in the M2c but not the M2a network in this study. IL-11 is a member of the IL-6-type cytokine family and may provide a link between inflammation and cancer [131, 132]. The role of IL-11 in macrophage activation has not been well studied yet.

Table 5 Boolean Functions and Description of IL-4 Stimulated M2a Activation Network. The * symbol denotes a non-validated link from public databases or literature
Table 6 Boolean Functions and Description of IL-10 Stimulated M2c Activation Network. The * symbol denotes a non-validated link from public databases or literature

M2a network has one input, IL-1α, which remains low consistently in the binarized profiles [71]. The M2c network has no input node. The shared output nodes for M2a and M2c include GM-CSF, IL-3, IL-6, IL-7, IL-12p70, IL-17, TNF-α, and OSM. IL-6 is an inflammatory cytokine and has demonstrated a consistently low profile in both M2 activation models [133, 134]. The low expression profile of IL-6 also indicates low expression profiles of Lymphotactin and IL-1a. TNF-α promotes the activation of inflammatory M1 cells and has a decreased value in M2 cells, which was also demonstrated by our M2 activation model [135, 136]. In addition, IL-6 and TNF-α, outputs of M2 macrophage, can reduce inflammation and kill the pathogen [134]. IL-12p70 is another common output that is elevated in M1 activation and remains as a logic low in both M2 models [137]. It has been reported that IL-17 induced the production of pro-inflammatory cytokines by human macrophage, but the quantitative analysis is still needed [138].

Validation of IL-4 induced M2a and IL-10 induced M2c activation models

Parent and child nodes, Boolean functions, attribution of each node, interaction strengths from linear programming algorithm, and published results confirming the links in M2a and M2c networks were also shown in Tables 5 and 6, respectively. The majority of the links can be verified based on the KEGG pathway or protein interaction database.

The established M2a Boolean network includes a link between MIP-2 (Chemokine (C-X-C) motif ligand 2-CXCL-2) and FGF-9 which cannot be validated by STRING or KEGG database. The regulation between MIP-2 and FGF-9 is the most complicated Boolean function in all three activation networks. However, it has been reported that FGF-9 enhanced M2 differentiation post-myocardial infarction accompanying with significantly elevated IL-10 secretion [139]. MIP-2 as an inflammatory protein is typically secreted by M1 macrophages. Another un-validated link in the M2a network is inhibition of IL-7 by MCP-5. Since MCP-5 is typically expressed in M1 and IL-7 in M2 [10, 11, 13, 136, 140], we consider these as meaningful Boolean regulations and keep them in our network model. In addition, all these unvalidated links are links to output nodes in the network, which will not affect the core regulation. Once the expression profiles in the experiments agree with the Boolean function, the links are kept in the established M2a activation network.

The only un-validated link in the M2c network is MIP-1β to FGF-9. Again, the experimental results illustrated the effectiveness of the Boolean regulation and the link was kept in the network model. FGF-9 demonstrated elevated expression levels in three macrophage activation networks, while LPS stimulated macrophages produced much more FGF-9 than M2a and M2c [22].

M2 activation models were further validated by temporal transitions of their Boolean network models. M2a and M2c activation networks were simplified by removing the output nodes and the single directional signal transduction nodes. The core network for M2a and M2c activation networks are represented by proteins in diamond nodes in Figs. 8 and 9. Simplified M2a network contains a 7-protein core network including FGF-9, MCP-1, MCP-5, MIP-1β, MIP-2, SCF, and VEGF-A. This core network has multiple feedback loops including MCP-1 attractor, MIP-2 attractor, MCP-1 and MIP-1β interaction loop, MCP-5 and VEGF-A interaction loop, and VEGF-A and SCF interaction loop. Besides, there are feedback loops among VEGF-A, MIP-2, FGF-9, SCF, and MCP-5 (Additional file 1: Figure S2).

Simplified M2c network contains a 5-protein core network including IL-11, IL-1α, SCF, MCP-5, and IFN-γ as shown in Fig. 9. This core network has MCP-5 and SCF self-inhibition loop, IL-1α, and IL-11 interaction loop, IL-1a and IFN-γ interaction loop, and feedback loops among IL-1α, IL-11, SCF, IFN-γ. The link between IFN-γ and IL-1α was obtained by removing MCP-3 since MCP-3 is just a signal transduction state in this branch. Separated core networks for M1, M2a, and M2c activation were also available in (Additional file 1: Figure S1, S2, and S3). All these interaction loops were confirmed by previous studies or the KEGG pathway as shown in Table 6.

Evolutionary maps of M2a and M2c were shown in Fig. 10 and Fig. 11, respectively. The evolution of M2a has a state transition map that contains 128 = 27 states while M2c has a state transition map for 32 = 25 states. The binarized M2a temporal profiles in our dataset showed the transition from \( {e}_{128}^3\to {e}_{128}^4\to {e}_{128}^{20}\to {e}_{128}^{84}\to {e}_{128}^{120} \) and these states were illustrated in red boxes in Fig. 10. The binary values of these 7 states were shown in Table 7. M2c network has binarized transitions from \( {e}_{32}^1\to {e}_{32}^5\to {e}_{32}^{13}\to {e}_{32}^{10}\to {e}_{32}^4\to {e}_{32}^{24}\to {e}_{32}^{32}\to {e}_{32}^{28} \) illustrated in red boxes in Fig. 11. The binary values of these 5 states were shown in Table 8. The computational binary state transition was compared with temporal experimental profiles and converged to experimental results.

Fig. 10
figure10

The evolution map of the core network for M2a activation showed all pathways predicted by the M2a activation Boolean network. The evolution of 128 (128 = 27) states for 7 nodes in the core network were illustrated from different initial conditions. The states in red boxes represented the evolution path based on the temporal cytokine profiles in this study

Fig. 11
figure11

The evolution map of the core network for M2c activation showed all pathways predicted by the M2c activation Boolean network. The evolution of 32 (32 = 25) states for 5 nodes in the core network were illustrated from different initial conditions. The states in red boxes represented the evolution path based on the temporal cytokine profiles in this study

Table 7 Binary Expressions of Proteins in the Core Network for IL-4 Stimulated Macrophage Activation. The binarization is shown in two’s complement representation
Table 8 Binary Expressions of Proteins in the Core Network for IL-10 Stimulated Macrophage Activation. The binarization is shown in two’s complement representation

Discussion

The most significant contribution of this study includes the following facts. We established and validated novel Boolean networks for 3 macrophage activation patterns using in vitro temporal cytokine profiles. Our results confirmed 1 stable M1 macrophage activation state and 2 stable M2 activation states, which agreed with reported in vitro classification on macrophage activation patterns. Validation of the models was conducted based on the coincidence of predicted cytokine expression levels, public databases, and reported experimental results. The literature review confirmed the majority of links in the Boolean networks for M1, M2a, and M2c network models as shown in Tables 3, 5, and 6, respectively. More importantly, our results demonstrated possible key factors IFN-γ, IL-1α, Lymphotactin, MCP-1, MCP-5, and SCF for M1 activation, FGF-9, MCP-1, MCP-5, MIP-1β, MIP-2, SCF, and VEGF-A for M2a activation IFN-γ, IL-1α, IL-11, MCP-5, and SCF for M2c activation. Further, expression profiles of these cytokines can serve as signatures for macrophage activation patterns since other protein expressions can be determined by proteins in the core network and input proteins of the network. These core signatures provide a possible mixed combination for macrophage activation markers instead of using only static individual markers. Additionally, cytokine markers might be more useful for in vivo study for easy and less costly measurements.

Robustness of Boolean networks is an important research aspect while only a few results have been reported [141,142,143]. These results have shown a large basin for the robustness and stability of a Boolean network. In this study, the simulated biggest core network for M1 macrophage activation has two stable states: \( {e}_{64}^7 \) and \( {e}_{64}^{58} \) as shown in Fig. 6. The state \( {e}_{64}^7 \) represents a state when MCP-5, Lymphotactin, IL-1α, and MCP-1 have elevated expression levels, SCF and IFN-γ have decreased expression levels. This is a single state attractor meaning the Boolean network starting from this initial state will stay here forever. All of the other 63 states of the network will evolve to the state of \( {e}_{64}^{58} \), which is a controversial state of \( {e}_{64}^7 \) (SCF and IFN-γ are elevated and other 4 proteins are down-regulated). Similarly, the core network for M2a has two stable states \( {e}_{128}^{88} \), and \( {e}_{128}^1 \) (Fig. 10). In the state of \( {e}_{128}^1, \) all 7 proteins have elevated expression levels. The core network of M2c activation has only one stable state \( {e}_{32}^{28} \) as shown in Fig. 11. All these have demonstrated the robustness of the established Boolean networks.

The STRING PPI database integrates several PPI databases and is the most widely adopted database for protein-protein interactions. We only use the STRING PPI database to confirm possible interactions among proteins in the established Boolean networks. Since the STRING PPI database does not provide directional information for protein-protein interactions, the directional information of our Boolean networks is generated by the linear programming algorithm computationally. The confirmation of directions in the network can be indirectly confirmed by 1) their positions in up- or downstream of pathways since cytokines can be either stimuli or production of a pathway; and 2) the temporal sequence of significantly elevated or inhibited expressions in biological experiments.

Conclusion

The proposed method can also be applied to large scale datasets to establish regulatory Boolean networks for other biological processes. Due to the limited sampling time points in most of the deposited experimental datasets, Boolean networks might be much more promising than traditional differential equation models, which require determination and estimation on a huge amount of parameters. In addition, with measurements at limited time points, deep learning might not be a good approach to attack specific biological processes at the current stage. It’s worth to mention, this method can be applied to a small set of key factors extracted from big data, and the following computational prediction and validation of the small networks are much easier than the big network. Currently, our software package (Matlab) can easily handle the evolution of a Boolean network with 11 nodes, leading to a total of 211 transition states using a normal workstation. For larger Boolean networks with more than 11 nodes, high-performance computing is more suitable.

We did notice some limitations of the proposed modeling method due to noisy/outlier functions. Outliers can still cause the linear programming to gather an extra connection, drop a connection or give a wrong sign to a connection. Additionally, a time delay can also cause problems in determining Boolean regulations depending on the amount of delay. For example, if protein A promotes protein B in a motif, up-regulation of protein A should give an instant high expression of protein B in the mathematical model. However, time delays and variations of temporal responses are very common in biological responses. One of our modeling assumptions is the Boolean network is fixed without changing nodes and changing Boolean functions in the response time period. Accordingly, Boolean functions are determined with respect to the whole experimental time period, not the single time point In this study. Therefore, small time-delays or variations of temporal response may affect the Boolean regulation at a short period of time, but won’t affect the Boolean regulation in the whole time period. However, if experimental results have a long delay on the up-regulation of protein B, incorrect Boolean functions may be selected. Finally, the model will be limited to the proteins that were measured in the experiment. Any other proteins not sampled were left out. Further studies will be conducted to integrate cross-platform data for a more complete macrophage activation model in the future.

Three Boolean network models were established to elucidate the regulation of M1, M2a, and M2c macrophage activations. Prediction of the Boolean network models agree with experimental results and validated the established models. Signatures of core cytokine profiles were determined, which provided a possible examination for macrophage activations based only on the cell productions.

Methods

Data Preparation: A total of 27 cytokines’ expression levels from macrophages isolated from C57BL/6 J male mice (Jackson Labs, Bar Harbor, ME) were collected at 7 time points, setpoint (0 h, 0 h), 0.5 h, 1 h, 3 h, 6 h, 12 h, 24 h from 4 biological groups (4 replicates/group) including 1 control group and 3 macrophage polarized groups: M1 induced by LPS, M2a induced by interleukin 4 (IL-4), and M2c induced by interleukin 10(IL-10) [22].

All proteins in the established networks are cytokines, which are small, nonstructural proteins, including interleukins, chemokines, interferons, and tumor necrosis factors [144]. A total of 34 cytokines have been reported in the literature and are secreted by different cell types including macrophages, T-cells, neutrophils, and fibroblasts. In the manuscript by Melton’s group, 24 proteins were reported. IL-1α, IL-6, and TNF-α were measured in their experiments and included in their dataset, but not reported because the authors thought these three cytokines were related to M2b activation [22]. All 27 cytokines secreted by macrophages were used to establish the Boolean network model in this study.

Due to the limited temporal measurements, a curve fitting algorithm was applied to enrich the dataset and this dataset served as inputs for a linear programming algorithm to find the most significant protein (node) and connections (links) for network structure. The same enriched data was also analyzed by an iterative K-means method to find the threshold for binarization. Thus, a network could be established where each node represents a protein and each link between the nodes represents an association of two proteins. Meanwhile, the expression level of each node was defined as a Boolean variable ‘0’ or ‘1’. A Karnaugh map was used to determine the Boolean regulatory functions among these nodes/variables to create a Boolean network model. Verification of the regulations was conducted based on published interactions from PPI databases and literature searches. The predicted evolution of each Boolean network model for a macrophage activation pattern was further validated with public databases and published results.

Curve fitting

The smooth spine algorithm in Matlab was applied to approximate continuous secretion profiles for 27 proteins. The algorithm was chosen due to its ability to construct an accurate curve with noisy data. By minimizing the squared error between the constructed data and experimental data at the sampling time points, optimized fitting curvature was obtained using the following equation [145]:

$$ p{\sum}_i{a}_i{\left({x}_i-s\left({t}_i\right)\right)}^2+\left(1-p\right)\int {\left(\frac{d^2s}{d{t}^2}\right)}^2 dt, $$
(1)

where p denotes the smoothing parameter controlling the level between the smoothness of the function and the error between the experimental data xi and the value of the fitting function s(ti) at chosen time points ti. Parameter ai represents the weight of the error, which is set to be the default value in MATLAB. Once the continuous fitting function s(ti) was gathered, we divided the 24-h time span into 50-time points with 40 being evenly sampled from the initial reading to the 9th hour and the last 10 samples coming from the 9th hour to the final 24. This was set because the majority of regulations occurred within the first few hours [19,20,21,22]. No negative value is allowed in the fitting algorithm to match with real biological situations.

Binarization

In a Boolean network model, protein expression is either up-regulated or down-regulated, therefore, binarization of the fitting function is needed. Any protein with similar expression levels during the 24 h time period is considered non-significant and will not be modeled. The iterative K-means method was used to cluster the expression levels for each protein and determine an activation threshold for each stimulus and experiment [146, 147]. The objective function for iterative K-means method using the absolute difference from the point to the center of the cluster was defined as:

$$ J={\sum}_{j=1}^k{\sum}_{i=1}^n\mid {x}_i-{C}_j\mid, $$
(2)

Where k is the cluster number, Cj is the centroid of the jth cluster, i represents each sample, and n is the total number of samples (50 in this study). We started with 8 clusters for iterative K-Means method and reduced to 2 clusters. The average of the mean distance from each cluster was used as the activation threshold for a specific protein. Any expression level below the activation threshold is defined as “0” and above the threshold is defined as “1” for binarization.

The threshold could be a value of 0 if there was no reading at any of the time points, or no significant change in protein expression levels. However, if a threshold of 0 will be considered in our modeling procedure, any non-zero low expression measurement will be considered active. To avoid such a conflict, a threshold of value “0” was changed to a low value claimed by the detection range of the biological experiment, such as 1 × 10− 6 in this study.

Linear programming algorithm

A network with connections among 27 proteins and each interaction strength was determined by linear programming method with an R package [148, 149]. In the network model, each protein and interaction among proteins is represented as a node with a link. Each node is assigned with an activity level xi ϵ R+ which is the expression level of each protein. Each link is assigned a weight, wj, i, form parent node j to the child node i. The weights wj, i ϵ R with wj, i > 0 representing promotion and wj,I < 0 for inhibition. Also, each protein has an activation threshold, δi, determined by the Binarization method above. The activity level of each node xi can be determined by its initial status \( {w}_i^0 \), the activity level of a parent node xj and contribution weight wj, i from the parent node as

$$ {x}_i={w}_i^0+{\sum}_{j\ne i}{w}_{j,i}{x}_j, $$
(3)

Considering biological networks are sparse, we minimize the sum of the absolute link weights \( \sum \limits_{i,j}\mid {w}_{j,i}\mid \), the bias term \( \sum {w}_i^0 \), and noise effects denoted by ζl. A random variable ζlϵ R0+, l = 1., …, L, is introduced to accounts for data variations in the experimental data. L represents the set of inactive proteins in the experimental data and coincides with the number of constraints that might violate the linear programming. The objective function was defined as:

$$ \mathit{\min}{\sum}_{i,j}\mid {w}_{j,i}\mid {x}_j+\frac{1}{Y}{\sum}_i{w}_i^0+\frac{1}{\lambda }{\sum}_l{\upzeta}_{\mathrm{l}} $$
(4)

subject to the constraints of

$$ \forall \left(i,m\right)\ where\ {x}_{i,m}\ge {\delta}_i:{w}_i^0+{\sum}_{j\ne i}{w}_{i,j}{x}_{j,m}\ge {\delta}_{i,m} $$
(5)
$$ \begin{aligned} \forall \left(i,m\right)\ where\ {x}_{i,m}&\ge {\delta}_i:{w}_i^0+{\sum}_{j\ne i}{w}_{i,j}{x}_{j,m}\\ &\ge 0 + {\upzeta}_{\mathrm{l}} \end{aligned}$$
(6)

where the constraints are defined for each experiment performed, 1,2,3,4, representing the 4 replicates in each experimental group, and for each time point from the enriched fitting data. Constraint (5) selected proteins up-regulating specific protein while constraint (6) looks for proteins down-regulating the protein. Effects of the data variation are minimized by the term \( \frac{1}{\lambda}\sum \limits_l{\upzeta}_{\mathrm{l}} \) in the objective function (4).

In the objective function, the production term, |wj, i|xj shows the contribution to the child node expression xi level from the parent node xj.Two parameters γ and λ are non-negative weights of the bias term and slack variables representing the penalties on the initial values and data variations in the objective function. For an ideal parameter, λ, 5-fold cross-validation was used with the smallest mean square error limiting the value of λ a maximum set of 2(xj,m), where σ is the variance of a variable. Parameter Υ in the objective function (4) represents the weights on the bias \( {w}_i^0 \). A bigger value of Υ (Υ = 100) in this study will put more penalty on the link of the proteins rather than bias. Minimizing this term helps to find the most significant connection to the child node, xi.

We further modified the linear programming algorithm to simplify the network structure. For each child node xi, a weighted contribution |wj, i|xj from a parent node (xj) is calculated. The maximum value of all weighted contribution from a parent node can be determined. For any node, if a value of |wj, i|xj is larger than 70% of the maximum value related to the node, the link from node j to node i is kept; In addition, orphan nodes and a node without significant changes in expression levels were not included in our model. The final output of the modified linear programming is a set of nodes and links with interaction strength wj, i. Positive wj, i means promotion while negative wj, i means inhibition.

Establish the Boolean logic function

Theoretically, to determine a Boolean regulation with 3 nodes, 8 (23) states representing the status of the 3 nodes are needed. It means at least 8-time points are needed and each time point represents a unique status of the network. Since we only have data at 7-time points, curving fitting with respect to time and resample the fitted data at different time points is reasonable to enrich the dataset. The 50 samples generated from the fitted data allows us to determine a Boolean function for up to 5 nodes ideally since 25 = 32 < 50. Considering the same status of the system may occur at different time points and most biological sub-networks have dual or triplet nodes, a maximum 5-node sub-network motif is expected with 50-time points in this study. It’s worth to mention, that the number of time points to be selected depends on how fast the dynamics evolve and the complexity of the network, is not the focus of this study and won’t be discussed in detail.

Based on assigning expression levels ‘0’ or ‘1’ to each node, a truth table of a motif including parents and child nodes can be established. Boolean regulation between parent and child node was determined by the Karnaugh map.

If there is a mismatch of promotion or inhibition between Boolean relation obtained from Karnaugh Map and linear programming or experimental profiles, then an extended 200 samples were re-examined (50 from each replicate and 4 replicates in each group) to determine the Boolean function instead of using the average. This can occur when one set of data has a large variation from the other 3 replicates, causing a skew in the average or when a time delay affects promotion and inhibition responses. The 200-sample-set was also used to determine if self-regulation could have been an option, which was mainly done on the results that always remained logic high, occurring once in the M1 cell and validated with experimental results.

Mathematical modeling and validation

The linear programming approach determined the nodes and possible connections in a network. Boolean regulations among nodes were determined by the Karnaugh map and represented by Boolean functions. With each node in the network as a Boolean variable X, the evolution of the status of a network can be written as:

$$ X\left(k+1\right)=M\ltimes X(k), $$
(7)

where k represents current time point and k +1 for next time point, X (k) and X (k + 1) represent current and future states of the network, M is the state transition matrix determined by Boolean functions and represents the semi-tensor product [150]. We and other research groups have established computational tools for Boolean network models that are available for free download [150,151,152]. .An example of determining a 3-node network is shown in Fig. 12. The 8 states, state transition matrix, and evolution of the 8 states for the 3-node network was computed with our own software package [151]. All codes associated with macrophage polarization can be found in our GitHub link (https://github.com/RicardoRamirez2020/Macrophage_Boolean_Network_LP).

Fig. 12
figure12

The Boolean functions and evolutionary map for a 3-node network were illustrated

Establishment of a core network

Dynamics of Boolean networks become more complicated with increasing numbers of nodes and Boolean functions, resulting in difficulties in understanding and analysis of the network. Given a network with N nodes, there exist 2N states of the network and the state transition matrix will be 2N by 2N dimension, indicating a huge computational cost. Simplification of a Boolean network for easy analysis has been conducted to reduce both the number of nodes and functions while still keeping the capability of predicting all states of the network. It’s worth to mention, output nodes of a network only have incoming signals and do not affect the regulations of the network. In addition, if multiple nodes are connected by single direction links, all intermediate nodes only transduce the signal with promotion (keep the signal) or inhibition (Not operation). These nodes with single direction connections can be simplified by keeping the starting and ending nodes and the corresponding Boolean function can be derived with basic logic operations. Therefore, by removing the output and intermediate transduction nodes, a core network can be established to represent a Boolean network with much simpler Boolean functions and fewer states for simulation. Meanwhile, with the determined states of a core network, states of removed nodes can be reconstructed.

Availability of data and materials

The datasets used and analysed during the current study are available from the corresponding author on reasonable request.

Abbreviations

FGF-9:

Fiboblast Factor 9

GM-CSF:

Granulocyte-Macrophage Colony-Stimulating Factor

IC:

Immune Complex

IFN-γ:

Interferon Gamma

IL-10:

Interleukin-10

IL-11:

Interleukin-11

IL-12p70:

Interleukin-12, p70

IL-17A:

Interleukin-17A

IL-1α:

Interleukin-1 Alpha

IL-2:

Interleukin-2

IL-3:

Interleukin-3

IL-4:

Interleukin-4

IL-6:

Interleukin-6

IL-7:

Interleukin-7

IP-10:

(CXCL10) C-X-C Motif Chemokine Ligand 10

KC-GRO:

(CXCL1) – C-X-C Motif Chemokine Ligand 1

LP:

Linear Programming

LPS:

Lipopolysaccharide

M1:

Classically Activated Macrophages

M2:

Alternatively Activated Macrophages

MCP-1:

(CCL2) Monocyte Chemoattractant Protein 1

MCP-3:

(CCL7) Monocyte Chemoattractant Protein 3

MCP-5:

(CCL12) Monocyte Chemoattractant Protein 5

MIP-1β:

(CCL4) C-C motif Chemokine Ligand 4

MIP-2:

(CXCL2) C-X-C motif Chemokine Ligand 2

mRNA:

Messenger Ribonucleic Acid

OSM:

Oncostatin M

PPI:

Protein to Protein Interactions

RANTES:

(CCL5) Regulated on Activation Normal T Cell Expressed and Secreted

SCF:

Stem Cell Factor

TGF- β:

Tumor Growth Factor Beta

TIMP-1:

Tissue Inhibitor of Metalloproteinases

TLR:

Toll-Like Receptor

TNF-α:

Tumour Necrosis Factor Alpha

VEGF-A:

Vascular Endothelial Growth Factor A

References

  1. 1.

    Arnold CE, Whyte CS, Gordon P, Barker RN, Rees AJ, Wilson HM. A critical role for suppressor of cytokine signalling 3 in promoting M1 macrophage activation and function in vitro and in vivo. Immunology. 2014;141(1):96–110.

  2. 2.

    Benjamin Emelia J, Virani Salim S, Callaway Clifton W, Chamberlain Alanna M, Chang Alexander R, Cheng S, Chiuve Stephanie E, Cushman M, Delling Francesca N, Deo R, et al. Heart disease and stroke statistics—2018 update: a report from the American Heart Association. Circulation. 2018;137(12):e67–e492.

  3. 3.

    Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68(1):7–30.

  4. 4.

    Atri C, Guerfali FZ, Laouini D. Role of human macrophage polarization in inflammation during infectious diseases. Int J Mol Sci. 2018;19(6):1801.

  5. 5.

    Parihar A, Eubank TD, Doseff AI. Monocytes and macrophages regulate immunity through dynamic networks of survival and cell death. J Innate Immun. 2010;2(3):204–15.

  6. 6.

    Heidt T, Courties G, Dutta P, Sager HB, Sebas M, Iwamoto Y, Sun Y, Da Silva N, Panizzi P, van der Lahn AM, et al. Differential contribution of monocytes to heart macrophages in steady-state and after myocardial infarction. Circ Res. 2014;115(2):284–95.

  7. 7.

    Tomlinson JE, Žygelytė E, Grenier JK, Edwards MG, Cheetham J. Temporal changes in macrophage phenotype after peripheral nerve injury. J Neuroinflammation. 2018;15(1):185.

  8. 8.

    Chen Y-J, Zhu H, Zhang N, Shen L, Wang R, Zhou J-S, Hu J-G, Lü H-Z. Temporal kinetics of macrophage polarization in the injured rat spinal cord. J Neurosci Res. 2015;93(10):1526–33.

  9. 9.

    Martinez FO, Gordon S. The M1 and M2 paradigm of macrophage activation: time for reassessment. F1000prime Rep. 2014;6:13.

  10. 10.

    Murray Peter J, Allen Judith E, Biswas Subhra K, Fisher Edward A, Gilroy Derek W, Goerdt S, Gordon S, Hamilton John A, Ivashkiv Lionel B, Lawrence T, et al. Macrophage activation and polarization: nomenclature and experimental guidelines. Immunity. 2014;41(1):14–20.

  11. 11.

    Murray PJ. Macrophage polarization. Annu Rev Physiol. 2017;79(1):541–66.

  12. 12.

    Ferrante CJ, Pinhal-Enfield G, Elson G, Cronstein BN, Hasko G, Outram S, Leibovich SJ. The adenosine-dependent angiogenic switch of macrophages to an M2-like phenotype is independent of interleukin-4 receptor alpha (IL-4Rα) signaling. Inflammation. 2013;36(4):921–31.

  13. 13.

    Porcheray F, Viaud S, Rimaniol AC, Léone C, Samah B, Dereuddre-Bosquet N, Dormont D, Gras G. Macrophage activation switching: an asset for the resolution of inflammation. Clin Exp Immunol. 2005;142(3):481–9.

  14. 14.

    Mosser DM, Edwards JP. Exploring the full spectrum of macrophage activation. Nat Rev Immunol. 2008;8(12):958–69.

  15. 15.

    Aurora AB, Porrello ER, Tan W, Mahmoud AI, Hill JA, Bassel-Duby R, Sadek HA, Olson EN. Macrophages are required for neonatal heart regeneration. J Clin Invest. 2014;124(3):1382–92.

  16. 16.

    Wang Y, Yang T, Ma Y, Halade GV, Zhang J, Lindsey ML, Jin Y-F. Mathematical modeling and stability analysis of macrophage activation in left ventricular remodeling post-myocardial infarction. BMC Genomics. 2012;13 Suppl 6(Suppl 6):S21.

  17. 17.

    Yabluchanskiy A, Ma Y, DeLeon-Pennell KY, Altara R, Halade GV, Voorhees AP, Nguyen NT, Jin Y-F, Winniford MD, Hall ME, et al. Myocardial infarction superimposed on aging: MMP-9 deletion promotes M2 macrophage polarization. J Gerontol A Biol Sci Med Sci. 2016;71(4):475–83.

  18. 18.

    Ma Y, Halade GV, Zhang J, Ramirez TA, Levin D, Voorhees A, Jin Y-F, Han H-C, Manicone AM, Lindsey ML. Matrix metalloproteinase-28 deletion exacerbates cardiac dysfunction and rupture after myocardial infarction in mice by inhibiting M2 macrophage activation. Circ Res. 2013;112(4):675–88.

  19. 19.

    Martin S, Zhang Z, Martino A, Faulon J-L. Boolean dynamics of genetic regulatory networks inferred from microarray time series data. Bioinformatics. 2007;23(7):866–74.

  20. 20.

    Rex J, Albrecht U, Ehlting C, Thomas M, Zanger UM, Sawodny O, Häussinger D, Ederer M, Feuer R, Bode JG. Model-based characterization of inflammatory gene expression patterns of activated macrophages. PLoS Comput Biol. 2016;12(7):e1005018.

  21. 21.

    Walter W, Alonso-Herranz L, Trappetti V, Crespo I, Ibberson M, Cedenilla M, Karaszewska A, Núñez V, Xenarios I, Arroyo AG. Deciphering the dynamic transcriptional and post-transcriptional networks of macrophages in the healthy heart and after myocardial injury. Cell Rep. 2018;23(2):622–36.

  22. 22.

    Melton DW, McManus LM, Gelfond JAL, Shireman PK. Temporal phenotypic features distinguish polarized macrophages in vitro. Autoimmunity. 2015;48(3):161–76.

  23. 23.

    Wang N, Liang H, Zen K. Molecular mechanisms that influence the macrophage m1-m2 polarization balance. Front Immunol. 2014;5:614.

  24. 24.

    Xaus J, Cardó M, Valledor AF, Soler C, Lloberas J, Celada A. Interferon γ induces the expression of p21waf-1 and arrests macrophage cell cycle, Preventing Induction of Apoptosis. Immunity. 1999;11(1):103–13.

  25. 25.

    Williamson M, Bingham B, Viau V. Central organization of androgen-sensitive pathways to the hypothalamic-pituitary-adrenal axis: implications for individual differences in responses to homeostatic threat and predisposition to disease. Prog Neuro-Psychopharmacol Biol Psychiatry. 2005;29(8):1239–48.

  26. 26.

    Lindner M, Thümmler K, Arthur A, Brunner S, Elliott C, McElroy D, Mohan H, Williams A, Edgar JM, Schuh C, et al. Fibroblast growth factor signalling in multiple sclerosis: inhibition of myelination and induction of pro-inflammatory environment by FGF9. Brain. 2015;138(7):1875–93.

  27. 27.

    Yoshihara S, Li Y, Xia J, Danzl N, Sykes M, Yang Y-G. Posttransplant Hemophagocytic Lymphohistiocytosis driven by myeloid cytokines and vicious cycles of T-cell and macrophage activation in humanized mice. Front Immunol. 2019;10:186.

  28. 28.

    Yamaguchi R, Yamamoto T, Sakamoto A, Ishimaru Y, Narahara S, Sugiuchi H, Yamaguchi Y. Roles of myeloperoxidase and GAPDH in interferon-gamma production of GM-CSF-dependent macrophages. Heliyon. 2016;2(2):e00080.

  29. 29.

    Halstead ES, Umstead TM, Davies ML, Kawasawa YI, Silveyra P, Howyrlak J, Yang L, Guo W, Hu S, Hewage EK, et al. GM-CSF overexpression after influenza a virus infection prevents mortality and moderates M1-like airway monocyte/macrophage polarization. Respir Res. 2018;19(1):3–3.

  30. 30.

    Yoshimura T, Imamichi T, Weiss JM, Sato M, Li L, Matsukawa A, Wang JM. Induction of monocyte Chemoattractant proteins in macrophages via the production of granulocyte/macrophage Colony-stimulating factor by breast Cancer cells. Front Immunol. 2016;7:2–2.

  31. 31.

    Hu X, Chakravarty SD, Ivashkiv LB. Regulation of IFN and TLR signaling during macrophage activation by opposing feedforward and feedback inhibition mechanisms. Immunol Rev. 2008;226:41–56.

  32. 32.

    Zdrenghea MT, Makrinioti H Fau - Muresan A, Muresan A Fau - Johnston SL, Johnston Sl Fau - Stanciu LA, Stanciu LA: The role of macrophage IL-10/innate IFN interplay during virus-induced asthma. (1099–1654 (Electronic)).

  33. 33.

    Lam D, Harris D, Qin Z. Inflammatory mediator profiling reveals immune properties of chemotactic gradients and macrophage mediator production inhibition during thioglycollate elicited peritoneal inflammation. Mediat Inflamm. 2013;2013.

  34. 34.

    Ferret-Bernard S, Pierre SX, Bach J-M. In vitro induction of inhibitory macrophage differentiation by granulocyte-macrophage colony-stimulating factor, stem cell factor and interferon-gamma from lineage phenotypes-negative c-kit-positive murine hematopoietic progenitor cells. Immunol Lett. 2004;91(2):221–7.

  35. 35.

    Kopydlowski KM, Salkowski Ca Fau - Cody MJ, Cody Mj Fau - van Rooijen N, van Rooijen N Fau - Major J, Major J Fau - Hamilton TA, Hamilton Ta Fau - Vogel SN, Vogel SN: Regulation of macrophage chemokine expression by lipopolysaccharide in vitro and in vivo. (0022–1767 (Print)).

  36. 36.

    Zdrenghea MT, Makrinioti H, Muresan A, Johnston SL, Stanciu LA. The role of macrophage IL-10/innate IFN interplay during virus-induced asthma. Rev Med Virol. 2015;25(1):33–49.

  37. 37.

    Recent Findings from Radboud University Provides New Insights into Rheumatoid Arthritis (Disease-Regulated Gene Therapy with Anti-Inflammatory Interleukin-10 Under the Control of the CXCL10 Promoter for the Treatment of Rheumatoid Arthritis). In: Gene Therapy Weekly. 2016: 49.

  38. 38.

    Vasquez RE, Soong L. CXCL10/Gamma Interferon-Inducible Protein 10-Mediated Protection against &lt;em&gt;Leishmania amazonensis&lt;/em&gt; Infection in Mice. Infect Immun. 2006;74(12):6769.

  39. 39.

    Mérida S, Sancho-Tello M, Almansa I, Desco C, Peris C, Moreno M-L, Villar VM, Navea A, Bosch-Morell F. Bevacizumab diminishes inflammation in an acute endotoxin-induced uveitis model. Front Pharmacol. 2018;9:649.

  40. 40.

    Xie WR, Deng H Fau - Li H, Li H Fau - Bowen TL, Bowen Tl Fau - Strong JA, Strong Ja Fau - Zhang JM, Zhang JM: Robust increase of cutaneous sensitivity, cytokine production and sympathetic sprouting in rats with localized inflammatory irritation of the spinal ganglia. (0306–4522 (Print)).

  41. 41.

    Albulescu R, Tanase C, Codrici E, Popescu DI, Cretoiu SM, Popescu LM. The secretome of myocardial telocytes modulates the activity of cardiac stem cells. J Cell Mol Med. 2015;19(8):1783–94.

  42. 42.

    Patel OV, Wilson Wb Fau - Qin Z, Qin Z: Production of LPS-induced inflammatory mediators in murine peritoneal macrophages: neocuproine as a broad inhibitor and ATP7A as a selective regulator. (1572–8773 (Electronic)).

  43. 43.

    Camelo S, Lajavardi L, Bochot A, Goldenberg B, Naud M-C, Brunel N, Lescure B, Klein C, Fattal E, Behar-Cohen F, et al. Protective effect of Intravitreal injection of vasoactive intestinal peptide–loaded liposomes on experimental autoimmune Uveoretinitis. J Ocul Pharmacol Ther. 2009;25(1):9–22.

  44. 44.

    Salmiheimo ANE, Mustonen HK, Vainionpää SAA, Shen Z, Kemppainen EAJ, Seppänen HE, Puolakkainen PA. Increasing the inflammatory competence of macrophages with IL-6 or with combination of IL-4 and LPS restrains the invasiveness of pancreatic Cancer cells. J Cancer. 2016;7(1):42–9.

  45. 45.

    Evans R, Fong M, Fuller J, Kamdar S, Meyerhardt J, Strassmann G. Tumor cell IL-6 gene expression is regulated by IL-1α/β and TNFα: proposed feedback mechanisms induced by the interaction of tumor cells and macrophages. J Leukoc Biol. 1992;52(4):463–8.

  46. 46.

    Rougier F, Cornu E, Praloran V, Denizot Y. Il-6 and IL-8 production by human bone marrow stromal cells. Cytokine. 1998;10(2):93–7.

  47. 47.

    Beq S, Rozlan S, Gautier D, Parker R, Mersseman V, Schilte C, Assouline B, Rancé I, Lavedan P, Morre M, et al. Injection of glycosylated recombinant simian IL-7 provokes rapid and massive T-cell homing in rhesus macaques. Blood. 2009;114(4):816.

  48. 48.

    Bonder CS, Finlay-Jones JJ, Hart PH. Interleukin-4 regulation of human monocyte and macrophage interleukin-10 and interleukin-12 production. Role of a functional interleukin-2 receptor γ-chain. Immunology. 1999;96(4):529–36.

  49. 49.

    de Waal MR, Abrams J, Bennett B, Figdor CG, de Vries JE. Interleukin 10(IL-10) inhibits cytokine synthesis by human monocytes: an autoregulatory role of IL-10 produced by monocytes. J Exp Med. 1991;174(5):1209.

  50. 50.

    Müller K, Bischof S, Sommer F, Lohoff M, Solbach W, Laskay T. Differential production of macrophage inflammatory protein 1γ (MIP-1γ), lymphotactin, and MIP-2 by CD4+ Th subsets polarized in vitro and in vivo. Infect Immun. 2003;71(11):6178–83.

  51. 51.

    Kasama T, Strieter RM, Lukacs NW, Lincoln PM, Burdick MD, Kunkel SL. Interleukin-10 expression and chemokine regulation during the evolution of murine type II collagen-induced arthritis. J Clin Investig. 1995;95(6):2868–76.

  52. 52.

    Lo U, Selvaraj V, Plane JM, Chechneva OV, Otsu K, Deng W. p38α (MAPK14) critically regulates the immunological response and the production of specific cytokines and chemokines in astrocytes. Sci Rep. 2014;4:7405.

  53. 53.

    Iwata H, Goettsch C, Sharma A, Ricchiuto P, Goh WWB, Halu A, Yamada I, Yoshida H, Hara T, Wei M, et al. PARP9 and PARP14 cross-regulate macrophage activation via STAT1 ADP-ribosylation. Nat Commun. 2016;7:12849.

  54. 54.

    Szymczak WA, Deepe GS. The CCL7-CCL2-CCR2 Axis regulates IL-4 production in lungs and fungal immunity. J Immunol. 2009;183(3):1964–74.

  55. 55.

    Wahl AF, Wallace PM. Oncostatin M in the anti-inflammatory response. Ann Rheum Dis. 2001;60(suppl 3):iii75.

  56. 56.

    Kashiwagi M, Hosoi J, Lai J-F, Brissette J, Ziegler SF, Morgan BA, Georgopoulos K. Direct control of regulatory T cells by keratinocytes. Nat Immunol. 2017;18(3):334–43.

  57. 57.

    Buchan G, Barrett K, Turner M, Chantry D, Maini RN, Feldmann M. Interleukin-1 and tumour necrosis factor mRNA expression in rheumatoid arthritis: prolonged production of IL-1 alpha. Clin Exp Immunol. 1988;73(3):449–55.

  58. 58.

    Herz J, Reitmeir R, Hagen SI, Reinboth BS, Guo Z, Zechariah A, ElAli A, Doeppner TR, Bacigaluppi M, Pluchino S, et al. Intracerebroventricularly delivered VEGF promotes contralesional corticorubral plasticity after focal cerebral ischemia via mechanisms involving anti-inflammatory actions. Neurobiol Dis. 2012;45(3):1077–85.

  59. 59.

    Mukherjee S, Chen L-Y, Papadimos TJ, Huang S, Zuraw BL, Pan ZK. Lipopolysaccharide-driven Th2 cytokine production in macrophages is regulated by both MyD88 and TRAM. J Biol Chem. 2009;284(43):29391–8.

  60. 60.

    Hu X, Chakravarty SD, Ivashkiv LB. Regulation of interferon and toll-like receptor signaling during macrophage activation by opposing feedforward and feedback inhibition mechanisms. Immunol Rev. 2008;226:41–56.

  61. 61.

    White AC, Lavine KJ, Ornitz DM. FGF9 and SHH regulate mesenchymal Vegfa expression and development of the pulmonary capillary network. Development. 2007;134(20):3743–52.

  62. 62.

    Teishima JUN, Yano S, Shoji K, Hayashi T, Goto K, Kitano H, Oka K, Nagamatsu H, Matsubara A. Accumulation of FGF9 in prostate Cancer correlates with epithelial-to-Mesenchymal transition and induction of VEGF-A expression. Anticancer Res. 2014;34(2):695–700.

  63. 63.

    Bhattacharya R, Ray Chaudhuri S, Roy SS. FGF9-induced ovarian cancer cell invasion involves VEGF-A/VEGFR2 augmentation by virtue of ETS1 upregulation and metabolic reprogramming. J Cell Biochem. 2018;119(10):8174–89.

  64. 64.

    Lam D, Harris D, Qin Z. Inflammatory mediator profiling reveals immune properties of chemotactic gradients and macrophage mediator production inhibition during Thioglycollate elicited peritoneal inflammation. Mediat Inflamm. 2013;2013:9.

  65. 65.

    Sierra-Filardi E, Nieto C, Domínguez-Soto Á, Barroso R, Sánchez-Mateos P, Puig-Kroger A, López-Bravo M, Joven J, Ardavín C, Rodríguez-Fernández JL, et al. CCL2 shapes macrophage polarization by GM-CSF and M-CSF: identification of CCL2/CCR2-dependent gene expression profile. J Immunol. 2014;192(8):3858.

  66. 66.

    Kim MS, Day CJ, Morrison NA. MCP-1 is induced by receptor activator of nuclear factor-κB ligand, promotes human osteoclast fusion, and rescues granulocyte macrophage Colony-stimulating factor suppression of osteoclast formation. J Biol Chem. 2005;280(16):16163–9.

  67. 67.

    Milks MW, Cripps JG, Lin H, Wang J, Robinson RT, Sargent JL, Whitfield ML, Gorham JD. The role of Ifng in alterations in liver gene expression in a mouse model of fulminant autoimmune hepatitis. Liver Int. 2009;29(9):1307–15.

  68. 68.

    Stoolman JS, Duncker PC, Huber AK, Giles DA, Washnock-Schmid JM, Soulika AM, Segal BM. An IFNγ/CXCL2 regulatory pathway determines lesion localization during EAE. J Neuroinflammation. 2018;15(1):208.

  69. 69.

    Reinders MEJ, Sho M, Izawa A, Wang P, Mukhopadhyay D, Koss KE, Geehan CS, Luster AD, Sayegh MH, Briscoe DM. Proinflammatory functions of vascular endothelial growth factor in alloimmunity. J Clin Investig. 2003;112(11):1655–65.

  70. 70.

    Gao N, Liu X, Wu J, Li J, Dong C, Wu X, Xiao X, Yu F-SX. CXCL10 suppression of hem- and lymph-angiogenesis in inflamed corneas through MMP13. Angiogenesis. 2017;20(4):505–18.

  71. 71.

    Dinarello CA. Interleukin-1 in the pathogenesis and treatment of inflammatory diseases. Blood. 2011;117(14):3720–32.

  72. 72.

    Jarmin DI, Nibbs RJB, Jamieson T, de Bono JS, Graham GJ. Granulocyte macrophage colony-stimulating factor and interleukin-3 regulate chemokine and chemokine receptor expression in bone marrow macrophages. Exp Hematol. 1999;27(12):1735–45.

  73. 73.

    Rull A, Beltrán-Debón R, Aragonès G, Rodríguez-Sanabria F, Alonso-Villaverde C, Camps J, Joven J. Expression of cytokine genes in the aorta is altered by the deficiency in MCP-1: effect of a high-fat, high-cholesterol diet. Cytokine. 2010;50(2):121–8.

  74. 74.

    Deshmane SL, Kremlev S, Amini S, Sawaya BE. Monocyte Chemoattractant Protein-1 (MCP-1): an overview. J Interf Cytokine Res. 2009;29(6):313–26.

  75. 75.

    Cotterell SEJ, Engwerda CR, Kaye PM. Leishmania donovani infection initiates T cell-independent chemokine responses, which are subsequently amplified in a T cell-dependent manner. Eur J Immunol. 1999;29(1):203–14.

  76. 76.

    Hogaboam CM, Lukacs NW, Chensue SW, Strieter RM, Kunkel SL. Monocyte Chemoattractant Protein-1 Synthesis by Murine Lung Fibroblasts Modulates CD4&lt;sup&gt;+&lt;/sup&gt; T Cell Activation. J Immunol. 1998;160(9):4606.

  77. 77.

    Ohta T, Sugiyama M, Hemmi H, Yamazaki C, Okura S, Sasaki I, Fukuda Y, Orimo T, Ishii KJ, Hoshino K, et al. Crucial roles of XCR1-expressing dendritic cells and the XCR1-XCL1 chemokine axis in intestinal immune homeostasis. Sci Rep. 2016;6:23505.

  78. 78.

    Alderson MR, Tough TW, Ziegler SF, Grabstein KH. Interleukin 7 induces cytokine secretion and tumoricidal activity by human peripheral blood monocytes. J Exp Med. 1991;173(4):923.

  79. 79.

    Wu W-K, Llewellyn OPC, Bates DO, Nicholson LB, Dick AD. IL-10 regulation of macrophage VEGF production is dependent on macrophage polarisation and hypoxia. Immunobiology. 2010;215(9):796–803.

  80. 80.

    Onishi RM, Gaffen SL. Interleukin-17 and its target genes: mechanisms of interleukin-17 function in disease. Immunology. 2010;129(3):311–21.

  81. 81.

    Croxford AL, Karbach S, Kurschus FC, Wörtge S, Nikolaev A, Yogev N, Klebow S, Schüler R, Reissig S, Piotrowski C, et al. IL-6 regulates neutrophil microabscess formation in IL-17A-driven Psoriasiform lesions. J Investig Dermatol. 2014;134(3):728–35.

  82. 82.

    Kelner GS, Kennedy J, Bacon KB, Kleyensteuber S, Largaespada DA, Jenkins NA, Copeland NG, Bazan JF, Moore KW, Schall TJ. Lymphotactin: a cytokine that represents a new class of chemokine. Science. 1994;266(5189):1395–9.

  83. 83.

    Son D-S, Roby KF. Interleukin-1α-induced chemokines in mouse Granulosa cells: impact on keratinocyte Chemoattractant chemokine, a CXC subfamily. Mol Endocrinol. 2006;20(11):2999–3013.

  84. 84.

    Vatakuti S, Schoonen WGEJ, Elferink MLG, Groothuis GMM, Olinga P. Acute toxicity of CCl4 but not of paracetamol induces a transcriptomic signature of fibrosis in precision-cut liver slices. Toxicol In Vitro. 2015;29(5):1012–20.

  85. 85.

    Zhang L, Zhao J, Kuboyama N, Abiko Y. Low-level laser irradiation treatment reduces CCL2 expression in rat rheumatoid synovia via a chemokine signaling pathway. Lasers Med Sci. 2011;26(5):707–17.

  86. 86.

    Kalinichenko VV, Bhattacharyya D, Zhou Y, Gusarova GA, Kim W, Shin B, Costa RH. Foxf1 +/− mice exhibit defective stellate cell activation and abnormal liver regeneration following CCl4 injury. Hepatology. 2003;37(1):107–17.

  87. 87.

    Ford J, Hughson A, Lim K, Bardina SV, Lu W, Charo IF, Lim JK, Fowell DJ. CCL7 is a negative regulator of cutaneous inflammation following Leishmania major infection. Front Immunol. 2019;9:3063.

  88. 88.

    Van Praag RME, Prins JM, Roos MTL, Schellekens PTA, Ten Berge IJM, Yong SL, Schuitemaker H, Eerenberg AJM, Jurriaans S, De Wolf F, et al. OKT3 and IL-2 treatment for purging of the latent HIV-1 reservoir in vivo results in selective Long-lasting CD4+ T cell depletion. J Clin Immunol. 2001;21(3):218–26.

  89. 89.

    Barnes PJ, Drazen JM, Rennard SI, Thomson NC. Asthma and COPD: basic mechanisms and clinical management: Elsevier; 2009.

  90. 90.

    Chiu B-C, Freeman CM, Stolberg VR, Komuniecki E, Lincoln PM, Kunkel SL, Chensue SW. Cytokine–chemokine networks in experimental mycobacterial and Schistosomal pulmonary granuloma formation. Am J Respir Cell Mol Biol. 2003;29(1):106–16.

  91. 91.

    Ramana CV, Gil MP, Han Y, Ransohoff RM, Schreiber RD, Stark GR. Stat1-independent regulation of gene expression in response to IFN-γ. Proc Natl Acad Sci U S A. 2001;98(12):6674–9.

  92. 92.

    Pype JL, Dupont LJ, Menten P, Coillie EV, Opdenakker G, Damme JV, Chung KF, Demedts MG, Verleden GM. Expression of monocyte chemotactic protein (MCP)-1, MCP-2, and MCP-3 by human airway smooth-muscle cells. Am J Respir Cell Mol Biol. 1999;21(4):528–36.

  93. 93.

    Shimizu Y, Irani AM, Brown EJ, Ashman LK, Schwartz LB. Human mast cells derived from fetal liver cells cultured with stem cell factor express a functional CD51/CD61 (alpha v beta 3) integrin. Blood. 1995;86(3):930.

  94. 94.

    Lin S-K, Chang H-H, Chen Y-J, Wang C-C, Galson DL, Hong C-Y, Kok S-H. Epigallocatechin-3-gallate diminishes CCL2 expression in human osteoblastic cells via up-regulation of phosphatidylinositol 3-kinase/Akt/Raf-1 interaction: a potential therapeutic benefit for arthritis. Arthritis Rheum. 2008;58(10):3145–56.

  95. 95.

    Migita K, Izumi Y, Torigoshi T, Satomura K, Izumi M, Nishino Y, Jiuchi Y, Nakamura M, Kozuru H, Nonaka F, et al. Inhibition of Janus kinase/signal transducer and activator of transcription (JAK/STAT) signalling pathway in rheumatoid synovial fibroblasts using small molecule compounds. Clin Exp Immunol. 2013;174(3):356–63.

  96. 96.

    Tsuboi I, Hirabayashi Y, Harada T, Hiramoto M, Kanno J, Inoue T, Aizawa S. Predominant regeneration of B-cell lineage, instead of myeloid lineage, of the bone marrow after 1 Gy whole-body irradiation in mice: role of differential cytokine expression between B-cell stimulation by IL10, Flt3 ligand and IL7 and myeloid suppression by GM-CSF and SCF. Radiat Res. 2008;170(1):15–22.

  97. 97.

    Preisler HD, Venugopla P, Sivaraman S, Larson R, Tricot G, Goldberg J, Miller K, Galvez A, Gregory SA, Adler S, et al. Selection of optimal remission consolidation therapy for individual patients with acute myelogenous leukemia. Exp Hematol. 2000;28(7):106–7.

  98. 98.

    Schook LB, Albrecht H, Gallay P, Jongeneel CV. Cytokine regulation of TNF-α mRNA and protein production by unprimed macrophages from C57BI/6 and NZW mice. J Leukoc Biol. 1994;56(4):514–20.

  99. 99.

    Nagoshi M, Sadanaga N, Joo H-G, Goedegebuure PS, Eberlein TJ. Tumor-specific cytokine release by donor T cells induces an effective host anti-tumor response through recruitment of host naive antigen presenting cells. Int J Cancer. 1999;80(2):308–14.

  100. 100.

    Bautz F, Rafii S, Kanz L, Möhle R. Expression and secretion of vascular endothelial growth factor-a by cytokine-stimulated hematopoietic progenitor cells: possible role in the hematopoietic microenvironment. Exp Hematol. 2000;28(6):700–6.

  101. 101.

    Dobrinski I. 750 Animal models to study germ line stem cells and spermatogenesis. J Anim Sci. 2017;95(suppl_4):364.

  102. 102.

    Larsson S, Löfdahl C-G, Linden M. IL-2 and IL-4 counteract budesonide inhibition of GM-CSF and IL-10, but not of IL-8, IL-12 or TNF-α production by human mononuclear blood cells. Br J Pharmacol. 1999;127(4):980–6.

  103. 103.

    Yamamoto Y, Klein TW, Friedman H. Involvement of mannose receptor in cytokine interleukin-1beta (IL-1beta), IL-6, and granulocyte-macrophage colony-stimulating factor responses, but not in chemokine macrophage inflammatory protein 1beta (MIP-1beta), MIP-2, and KC responses, caused by attachment of Candida albicans to macrophages. Infect Immun. 1997;65(3):1077–82.

  104. 104.

    Thomassen MJ, Raychaudhuri B, Bonfield T, Malur A, Abraham S, Barna B, Kavuru M. Elevated IL-10 Inhibits GM-CSF Synthesis in Pulmonary Alveolar Proteinosis. 2003;36.

  105. 105.

    Wilbers RHP, Westerhof LB, van Raaij DR, Smant G, Bakker J, Schots A. GM-CSF negatively regulates early IL-10 mediated responses. bioRxiv. 2017.

  106. 106.

    Shanmugam N, Reddy MA, Guha M, Natarajan R. High glucose-induced expression of Proinflammatory cytokine and chemokine genes in Monocytic cells. Diabetes. 2003;52(5):1256.

  107. 107.

    Orlofsky A, Wu Y, Prystowsky MB. DIVERGENT REGULATION OF THE MURINE CC CHEMOKINE C10 BY Th1AND Th2CYTOKINES. Cytokine. 2000;12(3):220–8.

  108. 108.

    Makita N, Hizukuri Y, Yamashiro K, Murakawa M, Hayashi Y. IL-10 enhances the phenotype of M2 macrophages induced by IL-4 and confers the ability to increase eosinophil migration. Int Immunol. 2015;27(3):131–41.

  109. 109.

    Gewiese-Rabsch J, Drucker C, Malchow S, Scheller J, Rose-John S. Role of IL-6 trans-signaling in CCl4 induced liver damage. Biochim Biophys Acta. 2010;1802(11):1054–61.

  110. 110.

    Natsume M, Tsuji H, Harada A, Akiyama M, Yano T, Ishikura H, Nakanishi I, Matsushima K. Kaneko S-i, Mukaida N: attenuated liver fibrosis and depressed serum albumin levels in carbon tetrachloride-treated IL-6-deficient mice. J Leukoc Biol. 1999;66(4):601–8.

  111. 111.

    White CA, Dimitriadis E, Sharkey AM, Stoikos CJ, Salamonsen LA. Interleukin 1 beta is induced by interleukin 11 during decidualization of human endometrial stromal cells, but is not released in a bioactive form. J Reprod Immunol. 2007;73(1):28–38.

  112. 112.

    Grigaitis P, Jonusiene V, Zitkute V, Dapkunas J, Dabkeviciene D, Sasnauskiene A. Exogenous interleukin-1α signaling negatively impacts acquired chemoresistance and alters cell adhesion molecule expression pattern in colorectal carcinoma cells HCT116. Cytokine. 2019;114:38–46.

  113. 113.

    Matsunaga K, Klein TW, Newton C, Friedman H, Yamamoto Y. &lt;em&gtLegionella pneumophila&lt;/em&gt; Suppresses Interleukin-12 Production by Macrophages. Infect Immun. 2001;69(3):1929.

  114. 114.

    Jacobsen SE, Veiby OP, Smeland EB. Cytotoxic lymphocyte maturation factor (interleukin 12) is a synergistic growth factor for hematopoietic stem cells. J Exp Med. 1993;178(2):413.

  115. 115.

    Ma J, Jung BG, Yi N, Samten B. Early secreted antigenic target of 6 kDa of mycobacterium tuberculosis stimulates macrophage Chemoattractant Protein-1 production by macrophages and its regulation by p38 mitogen-activated protein kinases and Interleukin-4. Scand J Immunol. 2016;84(1):39–48.

  116. 116.

    Elain G, Jeanneau K, Rutkowska A, Mir AK, Dev KK. The selective anti-IL17A monoclonal antibody secukinumab (AIN457) attenuates IL17A-induced levels of IL6 in human astrocytes. Glia. 2014;62(5):725–35.

  117. 117.

    Pérez CV, Pellizzari EH, Cigorraga SB, Galardo MN, Naito M, Lustig L, Jacobo PV. IL17A impairs blood–testis barrier integrity and induces testicular inflammation. Cell Tissue Res. 2014;358(3):885–98.

  118. 118.

    Chen K, Pociask DA, McAleer JP, Chan YR, Alcorn JF, Kreindler JL, Keyser MR, Shapiro SD, Houghton AM, Kolls JK, et al. IL-17RA is required for CCL2 expression, macrophage recruitment, and emphysema in response to cigarette smoke. PLoS One. 2011;6(5):e20333.

  119. 119.

    Feuser K, Thon K-P, Bischoff SC, Lorentz A. Human intestinal mast cells are a potent source of multiple chemokines. Cytokine. 2012;58(2):178–85.

  120. 120.

    Pattarini R, Smeyne RJ, Morgan JI. Temporal mRNA profiles of inflammatory mediators in the murine MPTP model of Parkinson’s disease. Neuroscience. 2007;145(2):654–68.

  121. 121.

    Yamamoto T. Pathogenic role of CCL2/MCP-1 in scleroderma. 2008;13.

  122. 122.

    Oliveira SHP, Lukacs NW. Stem cell factor and IgE-stimulated murine mast cells produce chemokines (CCL2, CCL17, CCL22) and express chemokine receptors. Inflamm Res. 2001;50(3):168–74.

  123. 123.

    Sanchez C, Deberg MA, Burton S, Devel P, Reginster J-YL, Henrotin YE. Differential regulation of chondrocyte metabolism by oncostatin M and interleukin-6. Osteoarthr Cartil. 2004;12(10):801–10.

  124. 124.

    Luhata PL, Moses Munkombwe N, Hatwiko H. Isolation and 1H-NMR identification of a tiliroside from Odontonema strictum (Acanthaceae). 2016;5.

  125. 125.

    Matsunaga T, Toba M, Teramoto T, Mizuya M, Aikawa K, Ohmori S. Formation of large vacuoles induced by cooperative effects of oncostatin M and dexamethasone in human fetal liver cells. Med Mol Morphol. 2008;41(1):53–8.

  126. 126.

    Kufe DW, Hait W, Holland JF, Frei E, Pollock RE. Holland-Frei Cancer Medicine 8, vol. 8: PMPH-USA; 2010. p. 686–709.

  127. 127.

    Wang J, Vodovotz Y, Fan L, Li Y, Liu Z, Namas R, Barclay D, Zamora R, Billiar TR, Wilson MA, et al. Injury-induced MRP8/MRP14 stimulates IP-10/CXCL10 in monocytes/macrophages. FASEB J. 2015;29(1):250–62.

  128. 128.

    Newby AC. Metalloproteinase expression in monocytes and macrophages and its relationship to atherosclerotic plaque instability. Arterioscler Thromb Vasc Biol. 2008;28(12):2108–14.

  129. 129.

    Nold M, Goede A, Eberhardt W, Pfeilschifter J, Mühl H. IL-18 initiates release of matrix metalloproteinase-9 from peripheral blood mononuclear cells without affecting tissue inhibitor of matrix metalloproteinases-1: suppression by TNFα blockage and modulation by IL-10. Naunyn Schmiedeberg's Arch Pharmacol. 2003;367(1):68–75.

  130. 130.

    Lacraz S, Nicod L, Galve-de Rochemonteix B, Baumberger C, Dayer JM, Welgus HG. Suppression of metalloproteinase biosynthesis in human alveolar macrophages by interleukin-4. J Clin Invest. 1992;90(2):382–8.

  131. 131.

    Xu DH, Zhu Z, Wakefield MR, Xiao H, Bai Q, Fang Y. The role of IL-11 in immunity and cancer. Cancer Lett. 2016;373(2):156–63.

  132. 132.

    Balic JJ, Garbers C, Rose-John S, Yu L, Jenkins BJ. Interleukin-11-driven gastric tumourigenesis is independent of trans-signalling. Cytokine. 2017;92:118–23.

  133. 133.

    Braune J, Weyer U, Hobusch C, Mauer J, Brüning JC, Bechmann I, Gericke M. IL-6 regulates M2 polarization and local proliferation of adipose tissue macrophages in obesity. J Immunol. 2017;198(7):2927–34.

  134. 134.

    Rőszer T. Understanding the mysterious M2 macrophage through activation markers and effector mechanisms. Mediat Inflamm. 2015;2015:Article ID 816460.

  135. 135.

    Zelová H, Hošek J. TNF-α signalling and inflammation: interactions between old acquaintances. Inflamm Res. 2013;62(7):641–51.

  136. 136.

    Opal SM, DePalo VA. Anti-inflammatory cytokines. Chest. 2000;117(4):1162–72.

  137. 137.

    Hsieh C, Macatonia S, Tripp C, Wolf S, O'Garra A, Murphy K. Development of TH1 CD4+ T cells through IL-12 produced by listeria-induced macrophages. Science. 1993;260(5107):547–9.

  138. 138.

    Jovanovic DV, Di Battista JA, Martel-Pelletier J, Jolicoeur FC, He Y, Zhang M, Mineau F, Pelletier J-P. IL-17 stimulates the production and expression of proinflammatory cytokines, IL-β and TNF-α, by human macrophages. J Immunol. 1998;160(7):3513–21.

  139. 139.

    Singla DK, Singla RD, Abdelli LS, Glass C. Fibroblast growth Factor-9 enhances M2 macrophage differentiation and attenuates adverse cardiac remodeling in the infarcted diabetic heart. PLoS One. 2015;10(3):e0120739.

  140. 140.

    Li R, Paul A, Ko KWS, Sheldon M, Rich BE, Terashima T, Dieker C, Cormier S, Li L, Nour EA, et al. Interleukin-7 induces recruitment of monocytes/macrophages to endothelium. Eur Heart J. 2012;33(24):3114–23.

  141. 141.

    Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci U S A. 2004;101(14):4781–6.

  142. 142.

    Willadsen K, Triesch J, Wiles J. Understanding robustness in random Boolean networks; 2008.

  143. 143.

    Bornholdt S. Boolean network models of cellular regulation: prospects and limitations. J R Soc Interface. 2008;5(Suppl 1):S85–94.

  144. 144.

    Dinarello CA. Historical insights into cytokines. Eur J Immunol. 2007;37(S1):S34–45.

  145. 145.

    MathWorks I: Smoothing Splines: Fit: Natwick : Math Works Inc., 1996.; 1996.

  146. 146.

    Mukherjee S, Hill SM: Network clustering: probing biological heterogeneity by sparse graphical models. Bioinformatics 2011, 27(7):994–1000.

  147. 147.

    Gupta JK, Singh S, Verma NK: Mtba: Matlab toolbox for biclustering analysis. In: IEEE Workshop on Computational Intelligence: Theories, Applications and Future Directions: 2013: IIT Kanpur India; 2013: 94–97.

  148. 148.

    Matos MRA, Knapp B, Kaderali L. lpNet: a linear programming approach to reconstruct signal transduction networks. Bioinformatics. 2015;31(19):3231–3.

  149. 149.

    Boyd S, Kim S-J, Vandenberghe L, Hassibi A. A tutorial on geometric programming. Optim Eng. 2007;8(1):67.

  150. 150.

    Cheng D, Qi H, Xue A. A survey on semi-tensor product of matrices. J Syst Sci Complex. 2007;20(2):304–22.

  151. 151.

    Liu R, Qian C, Liu S, Jin Y-F. State feedback control design for Boolean networks. BMC Syst Biol. 2016;10(3):70.

  152. 152.

    Rongjie L, Chunjiang Q, Yu-Fang J: Observability and sensor allocation for Boolean networks. In: 2017 American Control Conference (ACC): 24–26 May 2017 2017; 2017: 3880–3885.

Download references

Acknowledgments

The authors acknowledge the support by the National Science Foundation under grant no. 1826086, National Institutes of Health R01HL074236 and F30HL110743, Veterans Administration BX00186, and Valero Scholarship.

Funding

This work is supported in part by the National Science Foundation under grant no. 1826086, National Institutes of Health R01HL074236 and F30HL110743, Veterans Administration BX00186, and Valero Scholarship. The funders had no role in study design, data collection and analysis, interpretation of data, or preparation of the manuscript.

Author information

Affiliations

Authors

Contributions

RR, YFJ designed the research, RR, AMH, and JR conducted the research, DWM and PKS collected the data, RR, AMH, JR, CQ, and YFJ wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Yu-Fang Jin.

Ethics declarations

Ethics approval and consent to participate

All mouse studies and procedures complied with the National Institutes of Health Animal Care and Use Guidelines and were approved by the Institutional Animal Care and Use Committees of the University of Texas Health Science Center Protocol #01090E and the South Texas Veterans Health Care System Protocol #1107–002 (12025x).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ramirez, R., Herrera, A.M., Ramirez, J. et al. Deriving a Boolean dynamics to reveal macrophage activation with in vitro temporal cytokine expression profiles. BMC Bioinformatics 20, 725 (2019). https://doi.org/10.1186/s12859-019-3304-5

Download citation

Keywords

  • Macrophage polarization
  • Boolean networks
  • Cytokines
  • Inflammation