Combinatorial transcriptional codes specify macrophage polarization destinations

Macrophages are driven to form distinct functional phenotypes in response to different immunological stimuli, in a process widely referred to as macrophage polarization. Transcriptional regulators that guide macrophage polarization in response to a given trigger remain largely unknown. In this study, we interrogate the programmable landscape in macrophages to find regulatory panels that determine the precise polarization state that a macrophage is driven to. Towards this, we configure an integrative network analysis pipeline that utilizes macrophage transcriptomes in response to 28 different stimuli and reconstructs contextualized human gene regulatory networks, and identifies epicentres of perturbations in each case. We find that these contextualized regulatory networks form a spectrum of thirteen distinct clusters with M1 and M2 at the two ends. Using our computational pipeline, we identify combinatorial panels of epicentric transcription factors (TFs) for each polarization state. We demonstrate that a set of three TFs i.e., NFE2L2, BCL3 and CEBPB, is sufficient to change the polarization destination from M1 to M2. siRNA-mediated knockdown of the 3-TF set in THP1 derived M0 cells, despite exposure to an M1 stimulant, significantly attenuated the shift to M1 phenotype, and instead switched to an M2-like phenotype. Further, the siRNA-mediated knockdown of the 3-TF set rendered the macrophages hyper-susceptible to S. aureus infection, clearly showing the effect of the knockdown.

3 Background Macrophages, the key components of the innate immune system, are responsible for processes ranging from the maintenance of healthy tissues, phagocytic clearance of microbes, to tissue repair and remodelling (Aderem, 2003;Greenberg and Grinstein, 2002;Mosser et al., 2021;Murray and Wynn, 2011). They do most of these in a bid to restore normal physiology and act as a bridge between innate and adaptive immunity. In response to a range of stimuli, including pathogen insults, macrophages are directed to distinct functional phenotypes in a process widely referred to as macrophage polarization (Benoit et al., 2008;Murray, 2017). The most well-described phenotypes are the M1 and M2 polarization states, where M1 macrophages are pro-inflammatory in nature with microbicidal activity, while M2 is associated with an anti-inflammatory response. Stimulation of macrophages with lipopolysaccharides (LPS) or interferon-gamma (IFNγ) typically leads to an M1 state, whereas treatment with interleukins (IL4, IL13) polarize macrophages into the M2 state (Mills et al., 2000;Murray, 2017). M2 macrophages are further sub-categorised into M2a, M2b, M2c and M2d depending on the activation triggers and the resulting transcriptional changes (Mantovani et al., 2004;Rőszer, 2015;Yao et al., 2019). However, this model is said to be oversimplified as it fails to capture the complexity exhibited by different macrophage states in an in vivo scenario where the interaction between different cytokines, interleukins, and other effector molecules define the final polarized state (Atri et al., 2018;Martinez and Gordon, 2014) . A recent study on the transcriptomes of human monocyte-derived macrophages treated with a variety of physiological stimulants identified nine different polarization states, which better fit the description of a spectrum of states rather than a classical binary M1/M2 axis (Xue et al., 2014). The different states in the polarization spectrum present complex transcriptional profiles as seen from altered gene expression patterns in several hundred genes. The importance of transcriptional regulation in defining macrophage polarization has been brought out by several studies (Gerrick et al., 2018;Hörhold et al., 2020;Palma et al., 2018;Xue et al., 2014;Zhao et al., 2021). In vitro and ex-vivo studies have shown complete reversal of macrophage polarization states depending upon the molecules (pathogen-associated molecular patterns (PAMPs), host-derived damage-associated molecular patterns (DAMPs) or cytokines, etc) present in the microenvironment (Sica and Mantovani, 2012;Stout et al., 2005;Zhang and Mosser, 2008).
Several studies have shown that activation of transcription factors (TFs) such as signal transducers and activators of transcription (STAT1, STAT3, STAT4 and STAT6), nuclear transcription factor-κB (NF-kB), NF-E2-related factors (NRF1 and NRF2), kruppel-like factors (KLF4 and KLF6), 4 CCAAT-enhancer-binding proteins (C/EBP-α and C/EBP-β) and peroxisome proliferator-activated receptors (PPAR-α, PPAR-δ and PPAR-γ) are associated with macrophage polarization (Lawrence, 2011;Li et al., 2018;Tugal et al., 2013). But, the available information is insufficient to predict the polarization state that will result in a given condition. A possible reason is that a set of TFs together defines the state rather than individual TFs by themselves. While the triggers that lead to different polarization states are known, the specific set of TFs that nucleate the change or serve as control points to guide macrophage polarization is not known.
The precise polarization state attained by the macrophages in response to pathogenic triggers have a direct bearing on the outcome of the disease (Benoit et al., 2008;Koziel et al., 2009;Pidwill et al., 2021). In response to bacterial pathogens such as Staphylococcus aureus, macrophages induce the transcriptional activity of pro-inflammatory genes, activating the M1 program. The M1 phenotype is known to have an activated microbicidal machinery that is utilized in clearing such acute infections. However, prolonged or excessive inflammatory response by M1 macrophages can be pathological as it results in tissue damage and destruction. M2 macrophages produce antiinflammatory mediators and resolve the inflammation (Koziel et al., 2009;Pidwill et al., 2021). Not surprisingly, the skewing of macrophage polarization is observed in several cases of severe bacterial infections, sepsis, tumour progression in cancers, and a range of other inflammatory diseases (Ardura et al., 2019;Hamilton and Tak, 2009;Sica and Mantovani, 2012;Zheng et al., 2018).
Although our understanding of macrophage biology is rapidly increasing, the precise molecular controls and their associated processes involved in macrophage switching to different functional states are not well characterized. Identifying combinations of molecules that guide macrophage polarization will greatly facilitate a deeper understanding of how infections are contained, provide a basis to predict the precise phenotype that will be attained in different situations and serve as a foundation for precision diagnostics and treatment strategies focused on the macrophages.
In this work, we seek to interrogate the programmable landscape in macrophages to identify regulatory panels that govern the precise macrophage polarization state. Towards this, we configured a computational pipeline that uses an integrative graph-theoretical approach where differential transcriptomes of macrophages were mapped onto the knowledge-based genome-scale gene regulatory network (hGRN), computed connected sets of top-ranked perturbations, and identified regulatory barcodes defining different polarization states ultimately. We then experimentally evaluated specific TFs in the M1 regulatory barcode in a THP-1 derived macrophage cell line. We also probed the role of these specific TF sets in the S.aureus infected 5 human monocyte-derived macrophages and cell-line infection model. Together, we validated genes in the M1 barcode and demonstrated that knocking down specific TF sets in the barcode switches M1 to an M2 state. We further showed that such a state switch leads to hyper-susceptibility to S.aureus infection.

Materials and Methods
Transcriptomes of macrophage polarization states: We retrieved transcriptomes of macrophages stimulated with a diverse set of stimulants (GSE46903) from NCBI-GEO (Barrett et al., 2013;Xue et al., 2014). Briefly, the authors of this study purified monocytes from peripheral blood mononuclear cells (PBMCs) of healthy controls and subjected them to granulocyte-macrophage colony-stimulating factor (GM-CSF) or macrophage-colony stimulating factor (M-CSF) for 72 hours (hrs) to generate baseline macrophages (Mo). Then, they treated the Mo with one of the twenty-eight different stimulants for different time points ranging from 30 mins to 72 hrs and generated transcriptomes from each stimulant and each time point. The stimulants spanned across the broad spectrum of pro and anti-inflammatory triggers to generate a wide range of polarization states. The triggers and their end time-points (n=145) considered for this analysis are shown in Table S1. Further, to study the relevance of macrophage polarization in acute bacterial infections, transcriptomes of peripheral macrophages exposed to S. aureus were retrieved from NCBI-GEO (GSE13670) (Koziel et al., 2009).
Raw data from the respective transcriptome experiments were retrieved and pre-processed in a platform-specific manner using appropriate packages in R (Carvalho and Irizarry, 2010;Gautier et al., 2004). Probes below the detection limit in ≥ 80% of the samples in each dataset were filtered out, and multiple probes associated with the same gene were summarized by taking an average. We then performed differential analysis for the macrophages treated with stimulants compared to their untreated counterparts, using the limma package in R (Ritchie et al., 2015). Genes with fold change ≥ ±1.3 and a False Discovery Rate corrected q-value ≤0.05 computed using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995)were considered to be statistically significant differentially expressed genes (DEGs).

Reconstruction of human Gene Regulatory Network (hGRN):
We reconstructed a comprehensive knowledge-based human Gene Regulatory Network (hGRN), consisting of Transcription Factor (TF) to Target Gene (TG) and TF to TF interactions. To achieve this, we curated the experimentally determined regulatory interactions (TF-TG, TF-TF) associated with human TFs (Wingender et al., 2013) in this study. These interactions were obtained from the following resources: (a) literature curated resources such as the Human Transcriptional Regulation Interactions database (HTRIdb) (Bovolenta et al., 2012), Regulatory Network Repository (RegNetwork) (Liu et al., 2015), Transcriptional Regulatory Relationships Unravelled by Sentencebased Text-mining (TRRUST) (Han et al., 2015) , TRANSFAC resource from Harmonizome (Rouillard et al., 2016), (b) ChEA3 containing ChIP-seq determined interactions (Keenan et al., 2019) and (c) high confidence protein-protein binding interactions (TF-TF) were retrieved from human protein-protein interaction network -2 (hPPiN2) (Ravichandran et al., 2021). The resulting interactions from both categories were collated and made a non-redundant set. This resulted in a connected set of TF-TF and TF-TG interactions that comprise a 'master' human Gene Regulatory Network (hGRN) with 27, 702 nodes and 890, 991 interactions (Data file S1).

Model contexting and identification of epicentric TFset:
Starting from the hGRN, stimulantspecific GRNs were derived by weighting the master hGRN network with differential transcriptomes of stimulants (ie., stimulant versus baseline) and generated weighted networks for individual stimulants. A sensitive network mining approach, previously established in the laboratory was then applied to identify top-ranked activated and repressed sub-networks using Eq. 1, Eq. 2, and Eq. 3. This method of computing response networks which involves a knowledge-based network and a sensitive interrogation algorithm has been shown to outperform data-driven network inference methods in capturing biologically relevant processes and genes (Ravichandran and Chandra, 2019;Sambaturu et al., 2021). Briefly, the network mining algorithm works by computing minimum weight shortest paths, in which each path begins from a source node and ends with a sink node, identifying connected sets of edges that make up the least-cost paths. The shortest paths between all pairs of genes were computed using Dijkstra's algorithm using Zen library implemented in Python2.7. For a path of length n, the path cost was calculated as a summation of the edge weights EW ij of all edges forming the path, and normalized by the path length.
Subsequently, all shortest paths were sorted based on their path costs from lowest to highest, and the top 1% were considered to constitute the top-perturbed network (TPN). Nodes forming the TPN were evaluated for enrichment of DEGs using Fisher's exact test. Further, the epicentric nodes in each TPN were identified using the EpiTracer algorithm developed previously in the laboratory (Sambaturu et al., 2016). Epitracer estimates the influential state of a node by computing ripple centrality, which is the product of closeness centrality and outward reachability.
The TPN has two components, a top-activated network and a top-repressed network, generated by using appropriate node-weighting schemes. To generate the activated network, the weight of node i in a condition treated with Stimulant A (S A ) was computed as: Where FC was the fold change of gene i in a condition treated with stimulant A (S A ) with respect to the baseline S B (antilog values were used to compute fold changes). To generate the repressed network, the node weight of node i in SA was computed as: The edge weight EW ij in a given condition S A for an edge e comprised of nodes N i (S A ) and N j (SA) was calculated as Where N i (S A ) and N j (S A ) are the node weights of nodes i and j, respectively. Lower the edge weight, stronger is the interaction between the nodes.
These subnetworks were combined to obtain a top perturbed network (TPN) for each stimulant.
Nodes in the TPN were rank-ordered based on the ripple centrality score. Differentially expressed TFs with the ripple centrality score >0 were considered as influential transcription factors for a given TPN. The resultant TF-set defines the regulatory control for a given stimulant.

Clustering macrophage polarization states:
We performed Monte Carlo reference-based consensus clustering using the M3C package implemented in R (John et al., 2020) to identify the optimal number of macrophage polarization states. Spectral clustering with Pearson's correlation coefficient as a correlation metric for 200 iterations with 250 replicates was used to estimate the cluster structure that best fits the differential transcriptome of (a) entire transcriptome (n=12,164), and (b) TF-set1 (n=265). Further, Relative Cluster Stability Index (RCSI) and Monte Carlo adjusted p-values were computed using M3C to test the number and structure of the clusters. For the cluster visualization and comparison, Dendextend package in R (Galili, 2015) was used.

Network visualization and Enrichment analysis:
We have visualized networks using Cytoscape 3.2.0 in Allegro Spring-Electric layout (Shannon et al., 2003). We performed enrichment analysis for the genes of the regulatory cores of each cluster using Reactome pathway resource with default parameters (Fabregat et al., 2017) and the resultant hits with q-value ≤ 0.05 were considered to be significant (using a hypergeometric test).  Table S2.

Reagents, Cell line and
Transient transfection: PMA differentiated THP-1 cells were transiently transfected with 100nM targeted siRNA using low molecular weight polyethyleneimine (Sigma-Aldrich, 40872-7). Further, 36 hrs post transfection, cells were treated or infected as indicated for required time and processed for analysis.
In vitro CFU analysis: PMA differentiated THP-1 cells were infected with S. aureus Cowan I at an MOI of 10 for 2 hrs to facilitate internalization. Following 2 hours, cells were washed with complete DMEM (without antibiotics) and Gentamicin (100μg/ml) treatment was given for 2 hours to remove extracellular bacilli. Cells were lysed in 0.1% TritonX-100 containing 1X PBS and CFU estimation was done at 0hr (after gentamicin treatment) and 24hrs by plating appropriate dilution on LB plates.
Statistical analysis: Levels of significance for comparison between samples were determined by the Student's t-test distribution and one-way ANOVA followed by Tukey's multiple-comparisons.
The data in the graphs are expressed as the mean ± standard error for the values from at least 3 or more independent experiments and p values ≤ 0.05 were defined as significant. R version 3.6.3 was used for all statistical analyses.

Results
Subnetworks from contextualised hGRN identifies TF barcodes for macrophage polarization states: We first identified the 'programmable space' in macrophages in terms of the nature and range of transcriptional reprogramming that they achieve in response to different immunostimulants. The stimulants belong to broad classes of pro and anti-inflammatory cytokines such as Toll-like receptors (TLRs) ligands, immune complexes, interleukins, glucocorticoids or their selected combinations, comprehensively covering the macrophage polarization spectrum. We probed the programmable space by constructing and analyzing 28 unbiased 'response networks' of macrophages, corresponding to the treatment with 28 different immunostimulants. Towards this, we first constructed a knowledge-based human gene regulatory network (hGRN). Our hGRN covers regulatory elements from the whole genome and contains 27,702 genes and 890, 991 TF-TF and TF-TG directed interactions, which was rendered context-specific for each stimulant using their respective differential transcriptome data. We then subjected the 28 contextualized hGRNs to our network mining algorithm, which identified a connected set of genes that collectively vary the most in their expression levels in response to the given stimulant. This yielded 28 different subnetworks containing the highest ranked regulatory perturbations, with nodes ranging from 4,032 to 6,785 nodes (Data file S2). The resultant subnetworks containing connected sets of weighted and directed interactions capture cascades of regulatory events which are in essence the top-perturbed networks (TPNs) in each condition. Next, we configured a pipeline to identify the most influential transcription factors from the resultant TPNs. We used EpiTracer, an algorithm previously developed in the laboratory (Sambaturu et al., 2016) which identifies the highest ranked influential nodes in a given network using a ripple centrality measure. This step results in the identification of the genes that have the strongest 'causal' association with a given perturbation pattern. The word 'causal' here refers to the potential of the given gene for wielding an influence on the system, as estimated by the network analysis. Put together, we identified epicentric transcription factors (Ripple centrality score >0) which are differentially expressed in each of these TPNs. When pooled, these TFs from all stimulants, result in a set of 265 TFs (TF-set1 -Data file S3), defining the programmable space in relation to macrophage polarization.
Further, to probe the number of distinct polarization states, we carried out agnostic clustering based on the differential transcriptome of TF-set1, using a Monte Carlo reference-based consensus clustering (M3C) method. This resulted in thirteen different clusters. The clusters and the associated members are shown in (Table S3). The relative cluster stability index (RCSI) and Monte Carlo pvalues further validated that the cluster stability (Fig. 1A, B). Put together, it is clear that the TPNs group into 13 broad clusters, each representing a distinct polarization state (Fig. 1C). Further, the cluster member associations inferred from TF-set1 were found to be concordant with the number of clusters (Fig. S1A, B) and clustering pattern (Fig. S2, cophenetic coefficient = 0.68; p-value = 1.25e−51) inferred by considering the entire differential transcriptome (n=12,164). This clearly indicates that the TFs in TF-Set-1 are sufficient to characterize a given polarization state and that the subnetworks from the contextualized hGRNs define the macrophage polarization spectrum. similarly for all other states. These TF-sets serve as barcodes for each state where color gradient reflects the gene expression changes of each TF (as compared to an unstimulated state) (Fig. 1D).
Put together, these barcodes amount to a union of 223 TFs (TF-set2) that can be considered to form a pool of central regulatory factors for the polarization states studied here (TF-set2 -Data file S3).  We hypothesized that knock-down of NFE2L2, CEBPB, BCL3 in an M1 state should switch it to M2. To test this hypothesis, we designed an experiment (experimental design shown in Fig. 3A) where we take THP-1 derived macrophages (taken as Mo) and stimulate them with sLPS to obtain an M1 state. We first verified that sLPS treatment is leading to an M1 state by measuring known M1 markers (IL1B in Fig. 3B), which confirmed that they were indeed in an M1 state. We also measured the levels of M2 markers (IL10 in Fig.3B) and verified that there was no significant change. STAT4 is known to have an effect on the polarization state as its knockdown in M1 was shown to significantly attenuate the M1 state (Ishii et al., 2009) and we therefore used STAT4 as a positive control. We checked whether STAT4 was exhibiting the expected behaviour in our experimental setup and we found that was indeed the case (Fig. S3, Fig. 3C). We then performed siRNA mediated combined knock-down of NCB, the 3-TF set, to evaluate the M1 to M2 switch.
14 We verified that the siRNA mediated knockdown (referred to as siNCB) had worked by testing the gene expression values of the 3 NCB genes and finding them to be significantly reduced (Fig. S4,   Fig. 3D). Upon stimulation with sLPS in the siRNA treated cells, we observed that the known M1 markers were significantly reduced while the known M2 markers were significantly pronounced, clearly showing that the siRNA knockouts lead the cells to shift from M0 to M2 (Fig. 3D). The change was statistically significant when compared with the parent THP-1 cells without knockdown as well as with non-targeting siRNA controls. In summary, we observe that the knockdown of the 3gene set leads to a clear shift of the polarization state away from M1 towards M2, despite the sLPS trigger. THP1 monocytes were treated with PMA (20ng/ml) for 16h. NCB (NFE2L2, CEBPB , BCL3) specific siRNA was transfected. Post 8h, LPS (100ng/ml) treatment was given for 72h. Gene expression changes of M1 and M2 markers were quantified (Two-way ANOVA; P > 0.05 :ns, P < 0.05:*, P < 0.01:**, P < 0.001:*** ; N=3; NT: non-targeting siRNA)

Knockdown of NCB restricts the growth of S. aureus within macrophages
S. aureus is a facultative intracellular pathogen causing a broad spectrum of diseases such as minor skin infections to serious bloodstream infections and endocarditis. This pathogen is known to survive inside the host cell and escapes from professional phagocytic cells such as macrophages through various immune evasion strategies, especially through cytoprotective mechanisms by altering the host gene expression (Koziel et al., 2009;Pidwill et al., 2021).
We performed differential analysis of peripheral macrophages infected with S. aureus using the publicly available data (Koziel et al., 2009), and profiled the expression levels of TFs in M1 (C11) barcode and M1 markers. This analysis showed an upregulation of NCB and other TFs in the M1 barcode at the 8hr, 24 hr and 48hr time point (Fig. S5). We then tested the same in an in vitro S.aureus infection model using THP-1 macrophages and made a very similar observation, where the mRNA expression levels of the STAT1, NCB, M1, and M2 markers were quantified at 12hr, 24hr, 48hr and 72hr time points (Fig S6, S7). Specifically, the expression levels of NFE2L2, CEBPB, and BCL3 were significantly upregulated in the THP-1 macrophages upon S. aureus infection with respect to the uninfected control at 24 hrs (Fig. S7). This suggests that macrophages infected with S.aureus polarize to the M1 state within 24 hrs. Hence, this time point was used for further experiments.
Next, to evaluate the potential of NCB in regulating polarization of macrophages during S. aureus infection, we used the same in vitro infection setup but using the NCB knockdown cells. We observed a significant shift from M1 to M2 with the downregulation in the expression of M1 markers (CXCL2, IL1B, iNOS, SOCS3) (Fig. 4B, Fig. S8) and upregulation in the expression of M2 markers (TGFB and IL10) (Fig. 4B). Interestingly, we observed 1.6 fold significant increase in the burden of S.aureus within these macrophages (Fig. 4C). This clearly indicates that knockdown of NCB combination renders the macrophage more conducive for S. aureus to proliferate.

Discussion
Macrophages display a wide array of functions, from being killers of pathogens to healers of damaged tissue (Murray, 2017;Murray and Wynn, 2011;Wynn et al., 2013). Their versatility is attributed to their ability to respond to specific physiological or pathological cues, leading to their polarization into different states, each with a distinct functional phenotype. A range of molecular triggers that cause differentiation in macrophages is known (Hu et al., 2021;Jaguin et al., 2013;Xue et al., 2014). It is also well known that the differentiation spans a full spectrum of phenotypes with the classically activated M1 and the alternatively activated M2 states at two ends of the spectrum. The molecular profiles that define or cause these states, however, have not been sufficiently characterized. In particular, very little is known about the regulatory genes that control and guide the polarization in a specific direction. It is of great interest to identify TFs that are most influential in guiding the polarization. Given the importance of macrophages, several studies have investigated the variations in the transcriptomes in different diseases (Benoit et al., 2008;Mosser and Edwards, 2008;Ross et al., 2021;Schultze, 2016;Zheng et al., 2018) as well as upon exposure of macrophages to different exogenous and endogenous stimuli (Hu et al., 2021;Jaguin et al., 2013;Xue et al., 2014), which have resulted in a wealth of data that are publicly accessible.
These large datasets provide opportunities for addressing a variety of questions.
In this study, we algorithmically identify distinct polarization states, minimal regulatory signatures for each state as transcription factor barcodes and demonstrate that manipulation of the barcode leads to loss of switching, in an example case. Additionally, our work has indicated that arresting the switch to M1 leads to hypersusceptibility to S.aureus infection in cell lines, suggesting that M1 has a distinct role in limiting the extent of infection. Two aspects in our computational workflow are distinct from the previous methods, which facilitated the identification of barcodes characteristic of different macrophage polarization states. The first is the use of an unbiased knowledge-based gene regulatory network as compared to the more common correlation-based interaction networks. In the latter, the interactions are statistical associations, with no guarantee that the interaction is 'real' in the given biological system (Petri et al., 2015;Ravichandran and Chandra, 2019). We overcome this limitation by consolidating a comprehensively curated knowledge-based network where each interaction is based on experimental data reported in databases and in primary literature. We have previously carried out a systematic comparison and have shown that knowledge-based networks are superior to correlation-based networks in capturing biologically relevant genes and pathways (Ravichandran and Chandra, 2019). Our macrophage specific GRN (subset of hGRN) network contains 326, 281 interactions among 742 TFs and 11, 423 target genes, which is more comprehensive than the previously reported knowledge based macrophage GRNs in literature (Hörhold et al., 2020;Palma et al., 2018). The second is the network interrogation method. Our method of contextualizing the knowledge-based network for diverse stimulants using the corresponding transcriptomes allowed a quantitative comparison of variations over a common background network, yielding a top-perturbed regulatory subnetwork in each case. Our network interrogation method identifies the most influential 'epicentric' transcription factors, which lead to the derivation of quantitative barcodes of TF panels. To the best of our knowledge, such barcodes have not been reported earlier even for the well studied M1 and M2 states. We are hopeful our results of epicentric TFs and their molecular interaction trajectories will serve as key inputs for several investigations on the TFs individually as well as for specific TF panels. The response networks for each stimulant also serves as a rich resource of the top-ranked regulatory pathways in each case.
As a first step into the identification of barcodes, we used a pooled set of differentially expressed epicentric TFs (n=265) in different conditions that we identify from the networks and perform an 20 unbiased clustering of the macrophages. We used a rigorous clustering method, which produced a robust clustering pattern which clearly indicated that the 28 stimulants group into thirteen polarization states. A clustering pattern was described earlier for the same dataset based on the top 1000 most varying DEGs from each condition, which resulted in nine distinct states (Xue et al., 2014). The clustering pattern that we get here from the top-ranked epicentres in each, yielded a similar pattern, indicating that a subset of DEGs that occur in our topnets are sufficient to achieve the differentiation. In fact, our gene-regulatory models implied that the saturated (SA and PA) and unsaturated (LA, LiA, and, OA) fatty acids mediate different modes of resolution that were grouped together in the previous clustering pattern, and are now resolved into two sub-branches. Similarly, the effects of IFNγ and sLPS which were clubbed in the previous pattern are now resolved, consistent with the known regulatory differences (Hoeksema et al., 2015;Kang et al., 2019). The exact number of clusters is not the main message from this analysis. Instead, it is deriving a molecular basis for the differentiation of the functional states, of which M1 and M2 are two ends of the spectrum, while several other states are located in between the two in the polarization state spectrum, which can be described as a punctuated continuum.
Our results suggest that the barcodes could serve as the regulatory handles to attain different polarization states with different functional phenotypes. We test that hypothesis for the M1/M2 axis and indeed demonstrate that siRNA knockdown of the TFs in the barcode achieves polarization state switching from M1 to M2 polarization states. In the M1 barcode, we selected the TFs In conclusion, we report a systems level analysis of the gene regulatory networks to identify regulatory factors which define different macrophage activation states. Perturbation of these factors for reprogramming the macrophage subpopulations to desired polarization states can be explored for applications in controlling the inflammation gradient in different diseases such as sepsis, autoimmune conditions and chronic infections.
Author contributions: SR and NC conceptualized the study. KNB planned the validation experiments. SR carried out the data curation, computational methodology, analysis, and interpretation. BB carried out the experiments. NC and KNB acquired funding for the research.
NC supervised the whole project. SR and NC wrote the manuscript with inputs from KNB and BB.
All authors read and approved the final manuscript.

Declaration of interests: NC is a co-founder of qBiome Pvt. Ltd and HealthSeq Precision
Medicine Pvt Ltd, which had no role in this manuscript. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.    Supplementary tables Table S1. In vitro stimulant data. Table S2. The list of primers used for the qRT-PCT Table S3. The list of clusters, associated members, network details of regulatory cores and the number of TFs in the barcodes.

Supplementary data file
Data file S1: The human Gene Regulatory Network.

Data file S2:
Twenty eight different subnetworks capturing highest ranked regulatory perturbations.