To investigate further, we evaluated the standard deviation of Pearson r (0–1 h) for randomly chosen ORFs in steps of 10 up to 300 (with each selection repeated 30 times) from whole genome and measured the standard deviation (SD) of auto-correlation. We notice that the mean value of SD of auto-correlation decreases with number of ORFs selected (Figure 4E). Notably, at around 80 ORFs for wildtype (and 50 for single KOs) the auto-correlation transits to a flat trend (DKO, however, did not show any transition). The amplitude of SD of auto-correlation is large for less than 80 ORFs for wildtype (and around 50 for single KOs) and small otherwise (Figure 4E). These results reveal the coincidence of biphasic transition at around 50 ORFs. Therefore, we can consider acute response as acute-stochastic mode (high SD) and the weak response as collective mode (low SD) [27]. We also found biphasic response for auto-correlations between 1–4 h (Figure 4F). Investigating further the temporal response of acute-stochastic and collective modes we observed auto- and cross-correlation profiles of acute-stochastic mode is similar to immune-related genes whereas the average response of randomly selected 80 ORFs in collective mode is scaled to genome-wide response (Figure 4G–H, and Figure 3A), thus demonstrating that collective mode is scalable. This implies the collective mode corresponds to a global, a-specific adjustment of the expression network that can be appreciated even in the case of a sufficiently populated random sample of genes. The collective mode can be considered as a sort of ‘mean field’ encompassing the entire genome expression and needs a sample of at least 80 ORFs to be reliably detected. Emergence of regulatory signature from scattered expressions in all genotypes To understand the temporal progress of biphasic response, we investigated the changes of whole genome expression from early (0–1 h) to late (1–4 h) in each genotype. We plotted the expression change (Dx) of single ORFs (0–1 h for x-axis and 1–4 h for y-axis) (Figure 5A–D, left panels) and the average value of groups of ORFs (10, 50, 80, and 200 ORFs) (0–1 h for x-axis and 1–4 h for y-axis) sorted from highest to lowest expressions (Figure 5A–D, center panels for 200 ORFs and Figure S2). Remarkably, we observed the transition from large scatter in expression distribu- tions around the origin for single ORFs to smooth lines for group of 50 ORFs onwards for all genotypes. This is due to the fact that grouping of expression distribution for 50 ORFs or above forms Gaussian distribution and the average value of each group remarkably follows linear line (Figure S3A and see ‘Linear regressions analyses’ in Methods); the fact average values follow linear lines reveals the emergence of regulatory signature working at the level of groups of genes. Furthermore, the fluctuations on Gaussian distribution reduce as we increase the grouping size (Figure S3A–C). From these, we observe genes upregulated at early time points were downregulated at later time points and vice- versa, through switching in gene regulatory circuits. Biologically, upregulated ORFs in early signaling, when transcription rate is faster than mRNA decay, are downregulated by late signaling, when mRNA decay is prevalent. The vice versa occurs for downregulated ORFs. Notably, this switching behavior is also observed for DKO, thus reinforcing that DKO possess genome- wide LPS response through MyD88- and TRIF- independent manner. Deciphering gene regulatory mechanisms from emergent signature We observed earlier from auto and cross-correlations that the gene expression responses between genotypes are distinct from one Macrophage Genome-Wide Control another. To understand this in depths, we compared genome-wide expression changes between genotypes for 0–1 h, total of 6 combinations (Wildtype vs. MyD88 KO, Wildtype vs. TRIF KO, Wildtype vs. DKO, TRIF KO vs. MyD88 KO, TRIF KO vs. DKO and MyD88 KO vs. DKO) (Figure 6A–F). We plotted the expression change (Dx) of single ORF as well as taking the average value of each group of ORFs sorted from highest to lowest expressions. We observed i) large scatter in expression distributions around the origin for ORFs in the collective mode and ii) linear expression distribution of ORFs in the acute-stochastic mode (Figure 6A–F, left panels). These results further confirm biphasic response of LPS stimulation. For the collective mode, we found the distribution of averages from expression change (Dx) in each group follows transition from scatter to smooth lines for group of 50 ORFs onwards for all genotypes (Figure 6A–F, right panels & Figure S4), revealing the emergence of genome-wide control in wildtype and all KOs. This is due to the fact that when we group ORFs of 50 and above, the distribution of Dx becomes Gaussian and the mean values of distribution follow linearity in the same way as mentioned in previous section. Another observation is that DKO emergent signature is invariant from single KOs (Figure 6E–F, right panels). This further suggests that the contribution from DKO is independent from MyD88 and TRIF and derives from an unknown source (U). Thus, the emergent signature as smooth line can be used to decipher quantitative genome-wide differential control of transcriptional (Tr) and mRNA decay machineries (Dm) by considering individual roles of MyD88 (M), TRIF (T), their collective/synergistic roles (MT) and MyD88/TRIF-independent processes (U) for gene induction, i.e., wildtype = M+T+MT+U, MyD88 KO = T+U, TRIF KO = M+U and DKO = U. This linear superposition model (resembling analysis of variance scheme) cannot be applicable to analyze individual genes. However, for understanding averaging behaviors found in each genotype, it can help to shed light on overall control mechanisms of LPS stimulation. Wildtype upregulated ORFs in acute-stochastic mode (WTa+), are constituted of a single group of 80 ORFs whose average expression change in wildtype (M+T+MT+U) is 1.22, MyD88 KO (T+U) is 0.46, TRIF KO (M+U) is 0.80 and DKO (U) is 0. Hence, the relative contribution of each signaling to the activation of acute- stochastic mode ORFs is determined: M = 0.80, T = 0.46, MT = 20.04, U = 0. From these values, MyD88-dependant (M) and TRIF-dependant (T) pathways are important for transcription while MyD88 & TRIF-dependent synergized processes (MT) are insignificant (in contrast to collective mode, see below). Furthermore, unknown processes (U) are not indicated for acute-stochastic mode. Using DAVID analysis platform [28], the majority of biological processes of acute-stochastic mode is related to immune system, defense response, inflammatory response, etc., specifically activat- ed by M and T (Table 1, p,0.05). Next, we compared groups of upregulated ORFs in the collective mode of wildtype (WTc+) (x-axis) with MyD88 KO and DKO (y- axes) (Figure 6A and C, right panels). Unlike acute-stochastic mode, collective mode possesses average points in plot which follows flat profiles. This could be due to the equilibrium of transcriptional and mRNA decay machineries or shift in experimental or normalization process. To investigate this statistically, we selected grouping of equal number of random ORFs and assessed their average expression change in wildtype against all genotypes. Flat profiles were observed with M+T+MT+U = 0 (wildtype), T+U = 20.024 (MyD88 KO), M+U = 0.013 (TRIF KO) and U = 0.024 (DKO) (Figure 7A–C). Since random sampling displays genome-wide constant shift in MyD88 KO, TRIF KO and DKO, PLoS ONE | www.plosone.org 6 March 2009 | Volume 4 | Issue 3 | e4905