Bayesian modeling suggests that IL-12 (p40), IL-13 and MCP-1 drive murine cytokine networks in vivo
- Sarah L. Field†1,
- Tathagata Dasgupta†2,
- Michele Cummings1,
- Richard S. Savage3,
- Julius Adebayo2, 4,
- Hema McSara1,
- Jeremy Gunawardena2 and
- Nicolas M. Orsi1Email author
© Field et al. 2015
Received: 29 October 2015
Accepted: 31 October 2015
Published: 9 November 2015
Cytokine-hormone network deregulations underpin pathologies ranging from autoimmune disorders to cancer, but our understanding of these networks in physiological/pathophysiological states remains patchy. We employed Bayesian networks to analyze cytokine-hormone interactions in vivo using murine lactation as a dynamic, physiological model system.
Circulatory levels of estrogen, progesterone, prolactin and twenty-three cytokines were profiled in post partum mice with/without pups. The resultant networks were very robust and assembled about structural hubs, with evidence that interleukin (IL)-12 (p40), IL-13 and monocyte chemoattractant protein (MCP)-1 were the primary drivers of network behavior. Network structural conservation across physiological scenarios coupled with the successful empirical validation of our approach suggested that in silico network perturbations can predict in vivo qualitative responses. In silico perturbation of network components also captured biological features of cytokine interactions (antagonism, synergy, redundancy).
These findings highlight the potential of network-based approaches in identifying novel cytokine pharmacological targets and in predicting the effects of their exogenous manipulation in inflammatory/immune disorders.
KeywordsBayesian network Cytokine Hormone Lactation Mouse
Cytokines comprise an extensive array of extracellular protein mediators best known for their traditional immunoregulatory functions, although their multiple roles in the orchestration of an array of physiological and pathological processes such as cancer, autoimmunity and cardiovascular disease are now well-recognized [1–10]. Many former reductionist studies have attempted to ascribe given physiological effects to the actions of one or small numbers of cytokines, but this strategy has met with conceptual difficulties stemming from the fact that their actions can be paradoxical at different concentrations and differ according to the prevailing hormonal milieu [3, 11]. Furthermore, the mechanistic insight offered by such studies in terms of relative cytokine interactions is also limited given that physiological responses are seldom governed by any one cytokine but rather by the combined influences of many. In this respect, these mediators are believed to operate as part of highly complex networks, wherein they exhibit antagonism, synergy and functional redundancy [12–14]. It is therefore clear that gaining a fuller understanding of cytokine function in any biological or clinical scenario rests with clarifying their interactions at a network level rather than relying on the increasingly inadequate T helper cell type 1/2 (Th1/Th2) paradigm .
A Bayesian network can be inferred from experimental data through the correlations between experimentally-measured quantitative values of different nodes. Various machine learning techniques are used to undertake this inference process, which is often helped by the use of a prior-knowledge graph ‘seed’ incorporating well-recognized, literature-derived information which reduces the computational outlay required to learn networks from biological data. Such prior knowledge speeds up the search and avoids local minima, improving performance and yielding statistically more robust networks, as described in Djebbari and Quackenbush (2008) . Moreover, this bias does not limit the process to learn new interactions between the nodes. Accordingly, Bayesian networks are well adapted to noisy data, small sample sizes and, most importantly, a lack of detailed knowledge about how causal interactions are implemented at a biological level. Moreover, they also allow the effect of network perturbations to be explored in silico [18, 19, 22, 23]. Indeed, once the network structure is established, the impact of in silico perturbation of upstream nodes (i.e. by changing the values in the conditional probability bins for one or more nodes) can be tracked through the network structure to assess the changes in the conditional probabilities of the downstream nodes, thereby leading to specific, experimentally verifiable predictions. It is nonetheless possible to conceive more complex biological relationships which cannot be distilled into a Bayesian acyclic graph. In this respect, Bayesian networks are acyclic insofar as they do not allow for the portrayal of structural causal loops such as those which may arise from feedback relationships . However, they offer a sensible compromise between capturing complex causal relationships and computational feasibility, and have thus been used in a wide variety of scientific contexts [25–28].
While most Bayesian-based studies have focused on core biological processes using data from cell lines or tissue-based genome-wide expression microarrays, less attention has been paid to physiological or clinically relevant in vivo systemically-derived data [21, 29]. As such, Bayesian methodologies offer a valuable opportunity to define the nature of both known and novel mediator interrelationships in complex physiological processes such as immunity and inflammation. They are, however, correlative and their significance rests on the ability to draw causal inferences about the underlying biology from these relationships.
In order to validate the models learned through this methodology, both from a biological as much as from an analytical perspective, we developed two independent circulatory cytokine and hormone-based datasets: one drawn from lactating mice and the other from mice whose pups were removed at birth. With regard to the former network, murine lactation represents a unique model system to explore the dynamic physiological interactions between the hormonal environment and inflammation/immune function, with particular regard to the effects of prolactin (PRL), both a key driver of lactation and putative critical immunomodulator [30–33]. Since lactation PRL levels can conveniently be abrogated in vivo by removal of the suckling stimulus, this offers an excellent strategy to verify predictions made by in silico perturbation of the lactation network by comparing it with an independently-generated in vivo pup-free network with a physiologically-induced PRL abolition. These networks were further used to infer functionally significant nodes, such as ‘hubs’ and ‘drivers’, as well as characterizing the relationships between nodes, such as synergy and antagonism. To the best of our knowledge, there have been no previous attempts to employ Bayesian network analysis to unravel dynamic protein interactions based on physiological in vivo data in such a mechanistic and qualitatively validated manner (i.e. including in vivo confirmation of in silico perturbation-dependent predictions).
Most circulatory cytokines and hormones vary markedly throughout lactation
Cytokine and hormone concentrations throughout lactation
NC (n = 7)
Day 1 (n = 8)
Day 2 (n = 8)
Day 4 (n = 8)
Day 10 (n = 7)
Day 16 (n = 8)
Day 21 (n = 7)
Day 24 (n = 7)
0.52 ± 0.41a
3.40 ± 1.35ab
2.20 ± 0.75ab
7.46 ± 2.78bc
12.32 ± 2.42c
15.48 ± 2.16c
15.06 ± 2.85c
13.16 ± 1.81c
34.83 ± 11.89
50.12 ± 14.37
55.17 ± 6.45
57.21 ± 10.81
74.76 ± 6.21
51.07 ± 6.09
83.24 ± 25.51
41.59 ± 8.66
5.67 ± 1.41a
12.50 ± 2.43abc
10.87 ± 1.31ab
25.25 ± 3.68d
36.73 ± 6.09e
17.75 ± 2.77bcd
25.38 ± 4.81d
18.08 ± 3.49cd
1.50 ± 0.04a
1.90 ± 0.21ab
2.17 ± 0.51ac
2.77 ± 0.96ad
6.48 ± 3.33def
3.75 ± 0.84bdef
3.95 ± 1.13af
5.77 ± 2.36bcef
0.10 ± 0.03
0.21 ± 0.06
0.19 ± 0.03
0.24 ± 0.12
0.58 ± 0.13
0.22 ± 0.05
0.53 ± 0.25
0.36 ± 0.15
3.98 ± 0.69a
8.48 ± 1.50ab
10.32 ± 1.01bc
16.79 ± 2.54d
21.43 ± 2.24d
16.13 ± 2.48cd
17.84 ± 2.98d
13.79 ± 2.12cd
17.08 ± 7.19a
119.81 ± 25.30b
151.52 ± 56.55b
164.98 ± 37.82b
265.95 ± 40.50b
151.35 ± 20.92b
201.91 ± 51.93b
154.86 ± 41.57b
159.93 ± 24.69a
278.68 ± 28.05bc
236.61 ± 18.92acd
396.55 ± 46.52be
387.07 ± 28.31ef
394.43 ± 38.41be
296.10 ± 37.38bde
180.05 ± 20.59a
5.69 ± 2.78a
18.95 ± 6.92abc
15.56 ± 6.16ac
27.99 ± 5.74ad
109.05 ± 29.77d
31.03 ± 4.86cd
45.68 ± 12.24cd
105.96 ± 62.76bd
164.44 ± 16.75ab
137.59 ± 37.23abc
114.44 ± 8.51b
247.74 ± 31.16ac
213.84 ± 20.94ac
160.99 ± 27.04abc
232.08 ± 11.32cd
259.51 ± 28.02cd
42.34 ± 17.66a
316.15 ± 101.70bc
249.58 ± 46.91bd
682.30 ± 173.51be
1222.66 ± 326.6ef
786.09 ± 90.51ef
699.93 ± 84.99ce
912.65 ± 283.08cdf
117.81 ± 28.60a
324.07 ± 40.21b
261.46 ± 33.82ab
416.84 ± 67.09b
650.45 ± 108.10ab
334.77 ± 54.77ab
447.35 ± 86.24b
307.99 ± 69.89b
42.17 ± 11.33a
155.83 ± 54.39abc
243.46 ± 27.76bd
391.00 ± 62.99de
477.14 ± 51.54e
446.35 ± 26.19e
441.08 ± 66.12ce
424.34 ± 55.23e
115.26 ± 54.35
462.25 ± 158.06
388.22 ± 111.37
349.13 ± 55.09
517.32 ± 104.37
470.79 ± 42.03
446.18 ± 50.47
444.59 ± 33.41
1.28 ± 0.68a
25.06 ± 12.59bcd
11.39 ± 4.24ac
13.80 ± 2.29bc
40.38 ± 4.51d
20.45 ± 3.50bc
24.44 ± 5.15bcd
11.30 ± 2.23bcd
6.62 ± 2.31a
22.41 ± 2.29bc
20.15 ± 1.84b
30.83 ± 5.85cd
44.02 ± 4.48e
32.71 ± 2.76cd
37.76 ± 4.15de
27.38 ± 4.46cd
14.19 ± 3.16a
76.07 ± 14.31ab
51.22 ± 8.18a
150.66 ± 30.78cd
219.15 ± 44.86d
136.96 ± 14.43bc
160.27 ± 26.80cd
120.45 ± 35.03bc
9.95 ± 1.75a
19.05 ± 3.24ab
10.36 ± 1.51a
22.66 ± 2.22b
26.26 ± 3.91b
22.89 ± 1.65b
28.51 ± 3.82b
20.39 ± 3.47ab
109.21 ± 11.86a
277.14 ± 23.85bc
202.25 ± 18.47ab
311.17 ± 36.45c
440.78 ± 49.29d
295.60 ± 33.14c
327.16 ± 38.26c
259.00 ± 40.22bc
307.86 ± 44.49a
583.64 ± 48.62b
500.58 ± 38.00b
641.07 ± 66.25bc
813.92 ± 90.49c
501.03 ± 68.81b
672.14 ± 75.99bc
522.25 ± 72.06b
17.64 ± 8.90a
97.42 ± 19.81bc
85.99 ± 17.25c
152.07 ± 30.82cd
248.47 ± 22.27e
165.79 ± 11.07d
175.62 ± 25.75dbe
173.97 ± 28.40dbe
0.00 ± 0.00abc
2.21 ± 1.14ac
1.16 ± 0.96a
5.20 ± 2.82abc
4.79 ± 2.47abc
15.41 ± 3.52b
10.35 ± 2.82abc
10.73 ± 4.03bc
79.49 ± 41.61a
558.78 ± 154.41b
560.13 ± 188.57b
642.66 ± 169.61b
985.75 ± 148.06b
630.06 ± 43.70b
630.40 ± 125.34b
717.02 ± 173.20b
230.63 ± 40.27
274.75 ± 111.88
96.75 ± 15.47
159.81 ± 28.44
150.16 ± 27.68
123.11 ± 15.37
205.83 ± 23.55
129.20 ± 26.77
206.17 ± 46.07ab
106.89 ± 17.78a
241.39 ± 59.96ab
364.94 ± 58.13b
158.10 ± 52.16a
33.51 ± 5.28c
130.08 ± 47.18ac
243.94 ± 98.57ab
8.71 ± 4.72ab
78.88 ± 46.34acd
50.75 ± 19.27de
298.00 ± 79.05f
85.43 ± 17.26dfg
77.38 ± 43.77acdeg
1.14 ± 0.74b
7.00 ± 2.76bce
Cytokine and hormone profiles during lactation cluster into three distinct time-series
Development of Bayesian networks
Generation of the lactation prior network
In the search for a network graph, an initial prior network structure seed was used as a bias, and the final network was learned from the data. The space of possible networks was explored using the TabuSearch algorithm and, at each iteration, the score was evaluated by adding, removing or reversing individual edges. The introduction of a prior network seed initiated the learning stage from all 53 data samples for 26 nodes as soft bias wherein the original seed edges were judged by the score in the same way as any other network edge. For each node, all prior probabilities (as defined by the BDe scoring metric with Dirichlet distribution) exceeded zero, enabling the construction of Bayesian networks reflecting all possible nodal interactions. Initiating the parameters (conditional probabilities) for each node at the beginning was achieved using a uniform distribution; and these changed with the data. Details of the prior network (Additional file 2) revealed that this featured 21 nodes, 9 of which appeared as parents.
The lactation Bayesian network structure features six structural hubs
Sensitivity, specificity and F-score values for comparisons between lactation, pup-free and VBBSM networks
VBSSM vs Seeded lactation network
VBSSM vs Seeded pup-free network
VBSSM Lactation vs VBBSM Pup-free network
Exploration of networks
Perturbation of PRL in silico
Changes in conditional probability associated with perturbation of PRL and structural hubs
IL-13, MCP-1 and IL-12 (p40) are key driver nodes
As outlined above, IL-13 and MCP-1 were selected for perturbation by allocating them to a high concentration bin with a probability of 1 (based on their ~3 and 2-fold increases in high concentration conditional probabilities, respectively, following in silico PRL depletion) (Table 3). They were perturbed both individually and in combination. Increasing IL-13 concentration caused extensive network changes, including a shift in MCP-1 and GM-CSF to a higher concentration and IFN-γ to a moderately high concentration (Additional file 5 panel A). The effects extended as far downstream as the terminal node: KC became medium and MIP-1α, RANTES and TNF-α became more medium. Note that the term ‘medium’ is used herein as shorthand to refer to the intermediate, mid-range relative concentration ‘equal frequency’ bin based on the data discretization (i.e. the bin containing the third of samples falling in the middle of the range analyzed). Increasing MCP-1 resulted in similarly extensive changes: among its children, GM-CSF concentration increased markedly while MIP-1β and RANTES became medium (Additional file 5 panel B). Further downstream effects included MIP-1α concentration becoming high, KC moderately medium-high, and IFN-γ and TNF-α both more medium. PRL branch perturbation (combined PRL/IL-13 and PRL/MCP-1) also resulted in significant changes to downstream node conditional probabilities as far as the terminal node (Additional file 5 panel C). Similarly marked effects were noted for combined PRL/IL-13/MCP-1 perturbation, which resulted in marked changes in all downstream hub statuses and conditional probability values (Additional file 6).
As a major parent of the second network branch (i.e. the one without IL-3, E2 and eotaxin as first-line parents), IL-12 (p40) was chosen for perturbation. IL-12 (p40) perturbation in silico (by allocation to a high concentration bin) had a dramatic downstream impact: IL-1β, IL-12 (p70) and KC concentration increased, MIP-1α became moderately high, and MIP-1β, IFN-γ and TNF-α became medium (Additional file 7 panel A). Combined perturbations of both network branches (IL-12 (p40) and eotaxin) also affected common downstream nodes (Additional file 8, panels A and B).
MCP-1 can act synergistically with IL-13 and/or IL-12 (p40)
IL-3 and PRL are potentially antagonistic
Similar perturbations of the parent node IL-3 were also performed. Allocation of IL-3 to low, medium and high concentration bins affected IL-2, IL-13 and MCP-1 concentrations to some degree, although the most striking feature of these changes was that the greatest effects were noted when IL-3 concentration was medium, rather than high or low (i.e. its effects were not linearly related to concentration; Additional file 10). Perturbing IL-3 by allocating it to a high concentration bin resulted in an increase in IL-13 concentration, as did perturbing PRL to a low concentration. This suggested that, as parents, IL-3 and PRL may be antagonistic in terms of their effects on IL-13, a notion supported by the intermediate conditional probabilities for IL-13 concentration resulting from high levels of both these parents (Fig. 5).
MIP-1β exhibits a biphasic response to eotaxin
Allocating eotaxin to a high concentration bin had minor effects, causing little more than a shift in MIP-1β towards a higher concentration (Additional file 9 and Additional file 11 panel A). By contrast, allocating eotaxin to a low concentration bin had marked effects on its children: IL-9 concentration was reduced whereas MCP-1 and MIP-1β became more medium, an effect carried through downstream to GM-CSF (Additional file 11 panel B). Intriguingly, the shift in MIP-1β concentration was in the same direction, independent of whether eotaxin concentration was perturbed upwards or downwards, suggesting a concentration-dependent biphasic response.
In vivo validation
Cytokine profiles fall when lactation is not established
Comparison of Bayesian networks
Lactation and pup-free networks share striking similarities
As noted for the lactation network, the VBSSM variant for pup-free dams retained many of the core features of its Bayesian counterpart, with the exception of eotaxin which became a terminal node (Additional file 3, panel B). IL-12 (p40) and IL-10 were principal parents, with IFN-γ and MIP-1β displayed as structural hubs. As with the earlier case, the inference process was performed without bias of any prior knowledge of the interactions between the cytokines being studied. The associated F-core of 0.28 (as displayed in Table 2) is low, although this was expected given the relatively low number of edges and nodes observed (under very high confidence criteria) in comparison to the seeded network obtained in Fig. 4. The high true negative values and specificity (i.e. true negative rate) as well as the conservation of the major regulatory hubs again point to the relatively high topological similarity between the VBSSM and seeded networks (Table 2).
In silico perturbation can correctly predict in vivo responses
The most striking difference between the two physiological networks was that half of the nodes had moved to high concentration bins (a relative observation) in the pup-free scenario, echoing the node status distribution in the in silico perturbed network driven by high concentration IL-13 and MCP-1 (Fig. 7 and Additional file 5). Based on a comparison between predicted (Additional file 6) and monitored (Fig. 7) effects, PRL perturbation alone only categorized four nodes correctly/closely (P4, IL-1α, IL-2 and IL-10 to low concentration bins). However, combined PRL/IL-13/MCP-1 perturbation correctly classified the qualitative (i.e. quantitatively high/medium/low concentrations relative to the equal frequency bins) nodal status of P4, IL-1α, IL-10 (low concentration), RANTES, TNF-α (medium), IFN-γ (moderately high) and GM-CSF (high). Close categorization was noted for MIP-1β (predicted to be medium instead of moderately high), IL-6 (moderately low instead of medium), IL-17 (moderately low instead of low) and MIP-1α (moderately high instead of high) (Table 3 and Fig. 7). In addition to having MCP-1 as a parent, MIP-1α, MIP-1β and IFN-γ were also children of IL-12 (p40), and perturbing both of these parents by allocating them to a high concentration bin in silico led to their correct categorization in the pup-free network (Additional file 7 panel B and Fig. 7). Only the qualitative profile of G-CSF, KC, IL-5 and IL-9 in the pup-free network could not be predicted/validated by in silico perturbation of the lactation network (i.e. 79 % of nodal statuses were correctly predicted by the perturbations performed). Thus, the relative changes of levels (in both positive and negative directions ) as encoded in the probability tables broadly agreed with overall states of the nodes in the physiologically perturbed pup-free network model.
This study aimed to characterize physiological cytokine:cytokine and cytokine:hormone interactions in murine lactation as a model of inflammatory/immune mediator regulation using a Bayesian network-based approach. This method revealed a robust lactation network structure featuring two main branches organized around principal parents (IL-12 (p40), IL-3, E2 and eotaxin), structural hubs (IFN-γ, IL-13, MCP-1, MIP-1α, MIP-1β, and RANTES) and a terminal node (TNF-α). These pivotal roles resonate with the central role reported for these mediators in the control of cytokine networks: IFN-γ in atherogenesis , IL-12 (p40) in those governing T cell and macrophage responses , IL-13 in the pathophysiology of ulcerative colitis , IFN-γ, IL-12 (p40/p70) and MCP-1 in Erdheim-Chester disease  and eotaxin, IFN-γ, MIP-1α and MIP-1β in asthma . A second network was generated from data drawn from pup-free dams in which removal of the suckling stimulus physiologically prevented the rise in PRL levels characteristic of the onset of lactation. This second Bayesian network, which validated in vivo the predictions made in silico, had broadly similar topology (F score 0.861) and total network complexity and, most significantly, featured the same conserved structural hubs. These findings suggest that these hub nodes connected by high confidence edges may act as structural lynchpins around which both networks assemble, potentially playing a role in integrating diverse upstream cytokine-mediated signals to induce coordinated downstream responses. In silico network perturbation allowed a more detailed analysis of the relative contribution of various nodes to the control of the overall network behavior. The results also pointed to the presence of ‘driver’ nodes (nodes which, when perturbed, propagated maximum amounts of changes in terms of the number of downstream nodes affected). In this respect, one branch parent (IL-12 (p40)) and two structural hubs (IL-13 and MCP-1) were responsible for orchestrating the most significant changes within the network, thereby reinforcing the proposed concept of a signal integration role played by the latter two cytokines.
The original premise underlying the choice of lactation as a model system was that this physiological setting would feature marked changes in PRL, a hormone that has been proposed as a central regulator of cytokine networks . Surprisingly, PRL functioned neither as a parent nor as a structural hub, and its only connection to the rest of the network via P4 was through a weak confidence edge. However, the cluster analysis demonstrated that most cytokine profiles clustered separately and peaked later than PRL suggesting that these may be independently regulated in vivo, as reflected by the Bayesian analysis. Furthermore, when PRL was perturbed (reduced) in silico, there were, with the notable exception of P4, only minor qualitative shifts in conditional probability values in downstream mediators (including those for IL-13 - PRL’s connection to the rest of the lactation network - and MCP-1). These findings point to a limited role for PRL as an immunoregulator in murine lactation, the presence of extensive functional redundancy for its actions and/or, possibly, reflect the fact that its purported effects on cytokine profiles are largely drawn from the in vitro setting. Our findings are consistent with the observation that both PRL and PRL receptor knockout mice do not feature substantial immune dysfunction [46, 47]. The effects of PRL on MCP-1 are rather more unclear: in contrast to our findings, existing data implicate PRL as an inducer of MCP-1 in ovarian luteolysis and lactational bone resorption [48, 49], although the systemic/circulatory relationship between these mediators in response to the prolactinaemia of lactation remains unknown. Nevertheless, its effect on P4 supports PRL’s reported role in increasing serum P4 in rats .
In silico perturbation of downstream IL-13, by contrast, had marked effects on network response, suggesting that this cytokine had a threshold conditional probability at which downstream signaling was achieved which PRL alone as an IL-13 input function was unable to achieve. It is tempting to speculate that oxytocin, a lactation-related neurohypophyseal hormone which causes smooth muscle contraction during the let-down reflex, may represent an additional input given that the highest cytokine concentrations coincided with the period of maximal suckling/milk production [51–55]. This would be consistent with oxytocin’s known modulation of cytokine production and receptor-dependent activity in a range of settings [56, 57], and the ability of IL-13 to modulate oxytocin receptor expression .
A principal aspect of this study focused on the ability of in silico perturbation to predict cytokine levels in a different physiological scenario in vivo. This premise’s biological validation was achieved by comparing in silico lactation network behavior following perturbation induced by pup removal and the consequent abrogation of lactation in vivo. Strikingly, allocation of the driver nodes IL-13 and MCP-1 to high concentration bins in silico (akin to what occurred physiologically in the in vivo pup-free network) correctly (or closely) predicted the relative concentration status of 13 out of their 14 downstream cytokines in vivo. IL-13 and MCP-1 perturbation less accurately predicted MIP-1α and MIP-1β’s response to pup removal (based on their conditional probabilities). However, they were correctly categorized when their direct driver node parent (IL-12 (p40) and MCP-1) perturbations were used instead. These findings are consistent with the documented effects of MCP-1 induction of MIP-1α expression in murine aneurysm models and total IL-12 induction of MIP-1α in isolated human natural killer cells (albeit in the presence of IL-15) [59, 60]. By contrast, in silico perturbation failed to predict the cytokine concentration response of the peripheral terminal cytokine KC correctly. These findings may be accounted for by the possibility that despite KC being connected to its parents via high confidence edges in the lactation network, its regulation may be under the control of additional mediators and/or involve changes in cell-specific receptor expression not measured as part of this investigation.
In silico lactation network perturbation also pointed to a recognized feature of cytokine interactions: synergy (the joint action of two or more cytokines which, when acting in concert, potentiate each other’s effects). Previous studies have found IL-13 to be a selective inducer of MCP-1  as supported by the present data. Furthermore, we noted an additional and previously unreported synergistic relationship between these mediators in relation to IFN-γ concentration. When IL-13 and MCP-1 were perturbed together (i.e. both allocated to high concentration bins), this resulted in a greater increase in IFN-γ concentration than that observed when either was perturbed alone. Further, more striking evidence of synergy was noted following the combined perturbations of both network branches. IL-12 (p40) and MCP-1 perturbation independently increased the concentration of their direct child MIP-1α, but this effect was markedly greater when these parents were perturbed together. Analogous, though less pronounced synergistic interactions were also noted for combined IL-12 (p40) and IL-13 perturbation which, independently, are both known to increase IFN-γ concentration [62, 63]. To the best of our knowledge, these synergistic functional relationships have not previously been reported in the literature and will form the basis of our future investigations into how widespread/conserved these relationships occur across different physiological/pathological settings, body compartments (particularly in relation to immune privileged sites) and species.
In silico lactation network perturbation also revealed another functional property of inflammatory networks: antagonism (the opposing/modulatory action of one cytokine on another). In this regard, IL-3 and PRL had opposing actions on IL-13, with IL-3 seemingly increasing IL-13 concentrations (in line with its documented effects on cultured murine bone marrow cells ), such that when these parents were perturbed together (low PRL and high IL-3), the effect on IL-13 was greater than those induced by the independent perturbation of its parents (Fig. 5).
Another recognized property of cytokine networks is functional redundancy, wherein two independent mediators can fulfill the same role, thereby physiologically compensating for each other’s absence. In both in vivo scenarios, IL-4 was peripheral to the main Bayesian network despite tracking IL-13 levels, as established by the cluster analysis. This feature suggests functional redundancy, which is consistent with the fact that IL-4 and IL-13 are known to operate through the same receptor system [9, 14, 65]. It is worth noting that IL-6 and G-CSF were also orphaned in the pup-free dam network. It is unclear whether this reflects true functional redundancy rather than simply the smaller sample size used to construct this particular network. However, even in the lactation network where a larger overall sample size was used, the edges connecting IL-6 and G-CSF to the rest of the network were of low confidence. This highlights the caveat that the interpretation of Bayesian networks can benefit from additional a priori functional knowledge of the mediators investigated and that larger data sets - unsurprisingly - generate models which more accurately represent mediator inter-relationships, as demonstrated by the use of VBSSM models.
Another interesting phenomenon noted was the response of MIP-1β to changes induced by in silico eotaxin manipulation. Intriguingly, perturbing eotaxin to either high or low concentrations both resulted in MIP-1β concentration being higher. Concentration-dependent paradoxical effects are well recognized among cytokines and, in this instance, MIP-1β responses point to a concentration-dependent biphasic response to eotaxin. The implications of this observation are currently unclear but we speculate that this relationship allows scope for homeostatic control, possibly through the existence of physiological feedback loops undetectable using Bayesian methodologies. Alternatively, this may, as outlined for KC, reflect the influence of unknown mediators not measured as part of this investigation.
The broad conservation of cytokine relationships across physiological scenarios points to their integrated regulation and indicates that changes in a small number of driver nodes can potentially affect the concentration of multiple downstream mediators. Given their critical role in orchestrating cytokine networks, identifying driver nodes such as IL-12 (p40), IL-13 and MCP-1 may prove valuable in the exogenous manipulation of inflammatory networks. Major, more predictable network changes could thus be induced through the selective targeting of driver nodes (e.g. by using antibody-based interventions). Moreover, desired cytokine level modulation could be induced without causing major network disruptions in instances where terminal nodes are targeted instead. This would be consistent with current clinical practices such as the use of anti-TNF therapy in the management of a range of autoimmune disorders [66–68]. Interestingly, comparative trials of various agents in rheumatoid arthritis suggest that agents such as these (e.g. etanercept, certolizumab) are superior to those targeting other cytokines such as IL-1 (e.g. anakinra) , which reside further upstream in our models. If analogous networks could be constructed using clinical data, there would be scope for developing novel, targeted therapeutic interventions with more predictable immune-related side effects.
This Bayesian network-based analysis has proved valuable in clarifying the complex structure and causal murine systemic cytokine-hormone network relationships. However, these findings must be interpreted in the light of certain caveats. Firstly, the present Bayesian networks are necessarily incomplete by virtue of the fact that they do not contain the entire array of possible interacting mediators. Secondly, the interpretation of predicted qualitative changes in network behavior must be considered in a context-specific manner in order to glean physiologically meaningful inferences from the data (e.g. in relation to specific pathophysiologies or in a whole animal instead of cell-specific in vitro models which differ in functional receptor expression/cytokine production profiles). Thirdly, they do not represent all possible edges given that the methodology inherently precludes the use of structural loops . This accounts for why the present networks consistently feature TNF-α as the terminal node, despite studies indicating that, for example, TNF-α can induce MIP-1β expression in mice . Fourthly, these networks allow investigators to make valuable predictions about causality, although these may still require empirical verification through specific experiments in vivo and/or in vitro (e.g. by using specific exogenous cytokines, inhibitors/traps, antibodies, pharmacological agents, hormones, pathogens or physiological insults). Nevertheless, the present findings are, to the best of our knowledge, the earliest available data providing a detailed characterization and assessment of cytokine networks in a whole animal model. The networks generated were statistically robust and independently corroborated many established cytokine interrelationships described in the literature (e.g. between IL-13 and MCP-1, between IL-12 (p40) and IFN-γ) [61, 71].
The identification of synergy, antagonism, functional redundancy and concentration-dependent biphasic responses within these networks lifts this method of analysis from being purely descriptive to mechanistic. This suggests that Bayesian in vivo cytokine networks as shown herein describe real physiological changes in a non-biased fashion, in contrast to modeling endeavors performed on in vitro systems (which fall short of presenting a realistic physiological picture) or differential gene expression studies (which do not account for the post-transcriptional regulation of cytokine production) . Furthermore, the identification of conserved regulatory hubs points to the existence of a previously unknown core structure within these cytokine networks whose responses can be predicted with some accuracy. Whilst we believe that the exciting findings of this study are a significant first step towards improving our understanding of complex systemic inflammatory/immune networks, we remain mindful that many of the relationships described herein remain to be individually and more fully validated in the in vivo setting. We accept that multiple approaches may be needed in this context in order to build a comprehensive picture of multiple interactions, such as the use of (conditional) knock-out animal models and infusions of cytokine traps, antibodies, cytokines and their soluble/decoy receptors. To this end, our future work will focus on establishing whether the network structures that we have identified herein appear to be conserved across both a range of pathophysiological scenarios (e.g. cancer, autoimmune disorders, cardiovascular disease) as well as across species.
Figure 2 demonstrates the workflow used within this study. Initial biological analysis (Experiment 1, see below) investigated the cytokine/hormone profiles in lactating mice, from which Bayesian networks were developed and explored in silico. In vivo validation of networks was achieved by investigation of cytokine/hormone networks in non-lactating mice (Experiment 2).
Eight to ten-week old virgin CD1 female mice were group housed (10 per cage) with ad libitum access to water and Standard Beekay diet (B&K, Grimston, Aldborough, UK). The lighting cycle was 14 h:10 h light:dark, and humidity and temperature were maintained at 55-65 % and 21.5 ± 1 °C. Females were naturally pair-mated to 12–14 week old CD1 stud males of proven fertility following Whitten effect-induced estrus synchronization. Females were caged individually in late pregnancy to litter down and nurse their pups, then sacrificed throughout lactation on days 1 (<24 h of littering), 2, 4, 10, 16, 21 and 24 (n = 8, 8, 8, 7, 8, 7 and 7 animals, respectively i.e. 53 data points from which the lactation network was constructed; see later). Weaning occurred on day 21, when the independent pups were removed from their mothers. Samples were collected ±1 h half way through the lighting cycle to minimize the impact of circadian rhythms on any analytes measured. The number of pups per dam was adjusted to 8 by cross-fostering to standardize the suckling stimulus. Negative (baseline) controls were provided by naturally cycling virgin females of the same age and strain (n = 7).
An independent data set to test the predictive power (thereby providing biological validation) of the lactation Bayesian network was generated by preventing the establishment of lactation in dams whose entire litters were removed from them at birth (thus maintaining low PRL levels). These females were sacrificed on days 2 and 4 (time-matched) (n = 8 in both groups; these 16 data points were used to construct pup-free networks). Seventy-six mice were used in total.
The animals used in this study were sacrificed under Schedule 1 of the Animals (Scientific Procedures) Act, 1986 (UK). The use of different animals for each individual time point was required on both ethical and biological grounds given the severe physiological repercussions of collecting blood from lactating dams at such closely spaced time intervals. On one hand, this would have been inappropriate in causing unnecessary ‘pain, distress and lasting harm’ in the eyes of the UK legislative framework. On the other, significant, repeated blood loss would have upregulated both adrenocorticotropic activity and the production of haematopoietic endogenous colony stimulating factors, thereby affecting cytokine networks as a result of sampling rather than true changes in physiology. Moreover, dams subjected to such repeated stress would have been more prone to pup cannibalism, resulting in uneven suckling stimuli across females.
Sample collection and analysis
For the sake of text readability, cytokine and hormone acronyms are listed in the abbreviations section. Whole blood was collected by cardiac puncture as previously described , allowed to clot on ice and serum isolated by centrifugation at 5,000 rpm for 3 min. Serum was stored at −80 °C until analysis. The panel of cytokines chosen was based on the widest murine analytical array (with known immunoendocrine interactions) commercially available at the time of the study. Serum samples were analyzed for IL-1α, IL-1β, IL-2, IL-3, IL-4, IL-5, IL-6, IL-9, IL-10, IL-12 (p40), IL-12 (p70), IL-13, IL-17, eotaxin, G-CSF, GM-CSF, IFN-γ, KC, MCP-1, MIP-1α, MIP-1β, RANTES and TNF-α by multiplex immunoassay (Bio-Rad Laboratories, Hemel Hempstead, Hertfordshire, UK) on a Luminex-100 cytometer (Luminex Corporation, Austin, Texas), equipped with StarStation software (Applied Cytometry Systems, Dinnington, UK), as previously described . Samples were analyzed on one plate to avoid any potential batch or inter-plate variation. No missing values were identified, and samples falling below the level of detection of the assay were allocated a concentration of 0 pg/ml in relation to the blank in order to avoid skewing the data sets to an unrepresentative higher concentration. Hormones relevant to lactation with putative immunomodulatory effects were also selected for analysis: PRL concentrations were determined by homologous specific radioimmunoassay , while E2 and P4 were assayed by enzyme-linked immunosorbent assay according to the manufacturer’s instructions (Alpha Diagnostic, San Antonio, Texas).
Data presentation and analysis
Data were expressed as pg/ml (cytokines, E2) or ng/ml (PRL, P4) ± SEM. All data distributions were assessed for normality by Anderson-Darling tests. Basic analytical approaches were performed to highlight time course-related changes in analyte profiles and to better appreciate the data before applying machine learning approaches. These were based on subsequent Kruskall-Wallis/analysis of variance with post hoc Mann–Whitney-U/Fisher’s LSD tests, as appropriate. Pup removal data were similarly compared using t-tests or Mann–Whitney-U tests. Corrections for multiple comparisons were applied using the Benjamini and Hochberg False Discovery Rate method . Statistical analyses were performed using Minitab (Version 16) and ‘R’.
In order to identify significant changes in mean concentration of each cytokine over time throughout lactation, z-scores were computed for each of these, assuming a null hypothesis that the variable was constant at the weighted mean value, where the weight was (standard error)−2 i.e. the inverse square of the standard error (estimated for each cytokine/hormone at each time point using the available multiple measurements). Resultant z-score P values (0.05 threshold) were corrected using the False Discovery Rate. Time-series were also analyzed using Bayesian Hierarchical Clustering (BHC) in order to define the number of data clusters in a principled manner, wherein the BHC algorithm identified distinct groupings solely on the basis of input data . They were modeled as being drawn from one of a number of underlying curves and assigned to a specific cluster on strictly probabilistic grounds. Since the time-series were normalized to have zero mean and unit variance, the clustering analysis was sensitive only to their shape. A flexible, non-parametric regression Gaussian Process Model (herein inferring a nonparametric latent function over time) was fitted to the mean at each time point in order to model the trend underlying mean value changes over time for each mediator. Correlations between cytokine, steroid hormone and PRL profiles were determined using Pearson’s product–moment correlations as a basis for the heat map.
Bayesian network construction
In Bayesian network formalism, a network of interacting variables is represented as a graph in which the variables are nodes and their interactions are directed edges . The edge between two nodes, P1 and P2, is associated with a conditional probability table containing probability of the state of P2 given the state of P1. The approach only allows dependencies between a node and its immediate parents. Formally, a Bayesian network is defined to be a pair (G, ΘG) where G is a directed acyclic graph whose vertices are random variables Pi and ΘG is the conditional distribution for each variable given its parents: Pb(Pi | Pa(Pi)), where Pa(Pi) denotes the set of all parents of Pi in the graph. Conditional independence statements, encoded by the network structure, define the conditional probability distribution.
In order to establish a prior network containing proteins from the present analytical target set, a seed network learned from the biomedical literature, protein-protein interaction databases, or any combination thereof, was used. The information incorporated in the seed network was solely restricted to well-established mouse-based interactions, given that many cytokines exhibit species-specific differences in function owing to evolutionary adaptation of what is effectively a plastic system with inbuilt functional redundancy. MetaCore Inc. (GeneGo, a Thomson Reuters business, http://thomsonreuters.com/en/products-services/pharma-life-sciences/pharmaceuticalresearch/metacore.html) was principally used to generate the prior network along with thorough hand curation. The resultant seed network was then combined with the results obtained from ‘Predictionet’ (http://www.bioconductor.org/packages/devel/bioc/html/predictionet.html), a text mining web application which retrieves gene interactions reported in the literature by focusing on a core set of targets and the context within which the information pertaining to them is retrieved/based . Any conflicting edges in the prior network causing feedback cycles were removed since Bayesian networks inherently preclude the existence of structural loops: no node/child can be either its own ancestor or descendant (analogously, the prior of the network structure could not contain such loops either as it was incorporated as a multiplicative factor in the scoring metric). The network structure close to the prior network has higher probability: the parametric formula for this prior structure-related factor in the scoring metric is given in Heckerman et al. (1995) . Prior to performing the Bayesian network analysis, z-score normalizations were applied to the raw data in Matlab. A machine learning algorithm (implemented in the WEKA-based open-source package MeV ; http://www.cs.waikato.ac.nz/ml/weka/) was used to refine the seed network in conjunction with the experimental data derived from postpartum mice in order to predict a high-confidence network [21, 79]. We found that the literature mining and the protein:protein database (PPI) toolbox embedded within the MeV package proved largely ineffective for establishing a seed network.
Cytokine and hormone profiles were further discretized into categorical data for the BN analysis following z-score application in order to reduce computational expense, making these relative over time courses and thus a prerequisite for assessing network behavior, which is based on relative rather than absolute concentrations. These were assigned to three mutually exclusive equal width relative concentration bins i.e. rescaled/normalized data spreading between 0 and 1 were categorized into low, intermediate/neutral and high bands of different sizes but with the same frequency such that the number of samples allocated to each category was the same. Thus, in the learned networks, each protein had an underlying conditional probability table where the color of each node was determined by its own allied underlying histogram. The snapshot of a network thus represents relative (rather than absolute) concentrations of each protein in a given physiological setting (i.e. lactation or pup-free) following the independent categorization described above for each dataset.
Time-dependent autocorrelations may have been present across different temporal measurements in each series (since the data were drawn from a time course) such that there was a violation of the static BN learning assumption that training samples are independent. However, an assumption of independence between time points is less problematic in defining the 3 bin-based concentration states than the first-order Markov assumption in a typical dynamic state space model on continuous data. Moreover, the Spearman correlation coefficient between two nodes is significantly reduced when the system is randomly sampled over time with fewer time points. Hence, a static Bayesian network structure learning is expected to provide a more appropriate model for sparsely sampled data, as is often the case in a biological context. On the other hand, static modelling for a time series data is more meaningful under the assumption that the effects exist long enough to be visible across multiple time steps, such that causal relationships are sustained rather than one-off events. The Spearman correlation coefficients between two causally connected nodes remained high in both short and long sampling schemes. As regards data discretization, even in the presence of a periodic cycle, it would have been valid to apply this to concentration differences for the purposes of training a relationship model involving multiple variables. These data were then used to learn the Bayesian network. Both the network topology and edge-specific conditional probabilities were learned from the data, starting from the initial seed network. Only nodes with three parents were selected, as per convention in the field [21, 80]. Searching for a best possible network for a given set of moderate numbers of proteins or genes is computationally expensive. Among the various simplifying assumptions to tackle this problem, one is to reduce the number of parameters associated with the conditional probability tables for each non-root node by restricting the maximum number of parent nodes for a given child to be three - a reasonable assumption in a biological context. However, in the present study, the very high bootstrapping stringency resulted in more than three parents for some nodes in the final network. In this regard, standard non-parametric bootstrapping was applied (100 operations) to address potential over-fitting in the Bayesian analysis, wherein multiple data sets were created by re-sampling with replacement to estimate the confidence in the various network features, such as edges between the nodes learned . Several metrics and search algorithms were employed in combination in order to optimize the sensitivity of the results to the learning procedure scoring metric ; the Tabu Search algorithm was used to optimize the Bayes (BDe) score as the selected scoring metric.
Overall, 53 data points were used in total but none of the 81 cells filled by the prior probabilities were zero: this is inherent to the Bayesian metric with Dirichlet prior distribution (the Bayes and the special case BDe metric) structure. The scoring metric contains both the coefficients denoting the number of records in the dataset and the choices of priors on the counts coming from the Dirichlet distribution in order to fill all 81 cells. The assumptions behind the Bayesian measure of the goodness of a belief-network structure lead to strong constraints on the exponents so that all of them can be constructed . For a brief description of the Bayes and its special case the BDe metric, and the Tabu search algorithm used herein, please refer to the Weka manual: http://weka.sourceforge.net/manuals/weka.bn.pdf.
Bayesian networks represent a ‘snapshot’ at any given time; however, the networks generated herein reflected sustained effects or causal events seen across multiple time points. For Bayesian network generation, the interaction matrix between each node for each time point was learned and then used to initialize the inference process for the next time point. Therefore, the final Bayesian network visualized from the final time point reflected the dynamic process encoded by the previous time step. Overall, therefore, the structure of seeded networks was learned both from the data and prior knowledge network, which biased the search to a subspace. This way of optimization is an alternative to heuristic methods which avoids the likelihood of hitting local minima.
The robustness of the present approach was ensured by the bootstrapping experiment outlined above which provided good feature confidence measures and returned few low confidence edges, even after increasing the stringency of the bootstrap confidence to 0.9 (i.e. features occurring in ≥90 % of iterations). The network directed acyclic graph was then visualized using Cytoscape (http://www.cytoscape.org) .
VBSSM model generation
In order to further validate the robustness of this model, a VBSSM was built in Matlab . VBSSM implements an analytical approximation scheme to Bayesian state-space models and, unlike other related methods, does not take prior information (i.e. the seed network) into account for network reconstruction , such that networks were solely constructed from the experimental data. It should be noted that although time t + 1 involved 7–8 replicates independent from those of time t, these experiments used mice of the same strain/age with even litter sizes whose lactation started at approximately the same time under identical conditions. In this way, the sets of replicates across time points were assumed to be dynamically related so as to make VBSSMs meaningful for the purposes of internal, independent validation of the models learned using a seeded Bayesian network method.
VBSSM is a dynamic Bayesian network inference algorithm that uses linear Gaussian state-space models to help reverse-engineer interactions between proteins or transcriptional networks from time series data, thereby explicitly modelling the progression of a system over a multi-step time course. In state-space Bayesian models, the observed measurements depend on hidden states, which are assumed to evolve according to first order Markov chain. A variational approach to the above models presents a novel way to learn the structure and the optimal dimensions of the state space and, as such, the VBSSM algorithm provides distributions over the model parameters leading to an inference of its underlying structure. The majority of algorithms, including VBSSM, only allow edges across time steps and transition models that are stationary (linear time invariant) and first order Markovian. It is relatively efficient in terms of the size of the network search space and does not require sub computations, such as checking cycles within the network. However, the assumption of a stationary, first order Markov model is not appropriate for sparsely sampled datasets with a small number of time points each featuring a limited number of biological replicates as in the present data sets. For example, if node A inhibits B, the corresponding high negative Spearman correlation coefficient is significantly reduced when the system is randomly sampled over time with fewer time points. Hence, the traditional static Bayesian network structure learning is expected to provide a more appropriate model for sparsely sampled data (as is frequently the case in a biological context), which accounts for the use of static networks described above. Nevertheless, static modelling for a time series data is more meaningful under assumption that the effects exist long enough to be visible across multiple time steps and causal events are sustained causes, not one-shot triggers. The Spearman correlation coefficients between two causally connected nodes remain high in both short and long sampling schemes. Since the data are drawn from a time course, it invalidates the static Bayesian network learning assumption that training samples are independent. However, due to the nature of the sampling process, it is believed that a violation of this assumption is less problematic than the first-order Markov assumption.
In the present study, the VBSSM was constructed from 6 randomly chosen replicates from each time point (as there were 7–8 total replicates for each time point). The interaction matrix between nodes for each time step was learned and used to initialize the inference process for the next time step. In order to get the final static snapshot of the Bayesian network, the network was taken from the final time step, which takes into account the history of the dynamic process as encoded in the networks from previous time steps; recall, as outlined above, that VBSSMs are linearly time invariant.
The VBSSM algorithm allowed for the selection of a high significance level network interactions (measured in terms of z-score). In this regard, the z-score was used to calculate the normal cumulative distribution p, viz. that the probability that an observation from a normal distribution with mean zero and variance one will be less than z, giving a significance level of 1-p. Accordingly, a larger z-score means a lower significance level and a more stringent requirement level for the edges and nodes in a derived network . The networks obtained herein had a z-score of 2.33, which corresponded to a significance level of 0.01 and cumulative probability of 0.99.
Bayesian network perturbation and in vivo validation
Based on the results obtained in the first network, the systematic in silico perturbation of eotaxin, IL-3, IL-12 (p40), IL-13, MCP-1 and PRL nodes was performed by altering their conditional probability table values (thereby removing the uncertainties of concentration status) in order to determine their relative contribution towards the cytokine network structure . PRL (the key expected change between both experiments which would offer biological validation of in silico predictions) was perturbed first by allocating it to a low concentration bin (probability of 1) and the effects of this operation recorded as conditional probability tables which represented shifts in concentration (i.e. bin conditional probability) status of downstream nodes. This enabled us to determine its positive effects on immediate downstream hub nodes. Thereafter, as main first network downstream branch hubs from PRL, IL-13 and MCP-1 were perturbed (given a high concentration status), both individually and in combination. As a final step, IL-3 and IL-12 (p40) were perturbed alone and in combination with IL-13, MCP-1 and PRL in various combinations in order to determine their relative interactions and contributions to the network structure. The changes in conditional probabilities for downstream nodes elicited by these in silico perturbations were then compared to those obtained in vivo following pup-removal (which naturally abrogated lactation and thus prevented a rise in PRL).
Network structural comparisons
Where: specificity = TN/(TN + FP)
sensitivity = TP/(TP + FN)
As an additional measure of network similarity, the comparative complexity of these networks was also assessed as a representative measure given that this feature depends on both graphical structure and node categorization, and thus hinges on the number of parameters needed to define network probability distributions. The present Bayesian networks are categorical (i.e. discrete), with each node being assigned a value in a fixed set of categories. Total network complexity was thus measured as the sum of all of its node complexities, which were individually determined by the number of their respective parents and categories. A node with k parents with respective number of categories c1, c2, .., ck has complexity c1 x c2 x … x ck. Although in the original analysis the networks constructed were category 3 with 3 bins for evaluating each node’s conditional probability, category 2 was chosen herein during the simulation process for simplicity. This change in category did not affect any of results obtained because only structural features, in particular topological similarity between the different networks, were being investigated.
It is worth noting that the discretized concentration distributions and binning were independent for the lactation and pup-free data sets. This was deliberate given that, firstly, relative concentrations define a system’s behavior rather than absolute values , making it more meaningful to compare relative nodal statuses across models, especially if covering different time courses (and thus tailored to different biological endpoints as herein). Secondly, for unbiased model validation, it is crucial that the validation data set should not inform model development based on the training data set (a paradigm which would be have been violated had the categories been defined on the basis of pooled data).
Bayesian hierarchical clustering
- E2 :
granulocyte-colony stimulating factor
granulocyte macrophage-colony stimulating factor
monocyte chemoattractant protein
macrophage inflammatory protein
- P4 :
regulated upon activation, normal T cell expressed and secreted
standard error of the mean
tumor necrosis factor
Variational Bayesian state space model
Invaluable advice was provided by Dr Ian Overton, Dr Rebecca Ward, Dr Megha Padi, Dr Benjamin Haibe-Kains, Dr Leon Peshkin, Mr Raktim Sinha and Mr Taylor Southwick. The authors would particularly like to thank Dr Amira Djebbari and Professor John Quackenbush for invaluable discussions and support during manuscript preparation and revision. The authors are also grateful to Dr Albert Parlow (National Hormone & Pituitary Program) for providing the PRL assay service, and to Ms Katherine Martin and Ms Uma Ekbote for excellent technical assistance.
Open AccessThis 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.
- Crivellato E. The role of angiogenic growth factors in organogenesis. Int J Dev Biol. 2011;55(4–5):365–75.View ArticlePubMedGoogle Scholar
- Kopf M, Bachmann MF, Marsland BJ, Gouwy M, Struyf S, Proost P, et al. Averting inflammation by targeting the cytokine environment. Nat Rev Drug Discov. 2010;9(9):703–18.
- Lopez-Pedrera C, Perez-Sanchez C, Ramos-Casals M, Santos-Gonzalez M, Rodriguez-Ariza A, Cuadrado MJ. Cardiovascular risk in systemic autoimmune diseases: epigenetic mechanisms of immune regulatory functions. Clin Dev Immunol. 2012;2012:9746–8.View ArticleGoogle Scholar
- Ma S, Ma CC. Recent development in pleiotropic effects of statins on cardiovascular disease through regulation of transforming growth factor-beta superfamily. Cytokine Growth Factor Rev. 2011;22(3):167–75.PubMedGoogle Scholar
- Nowsheen S, Aziz K, Panayiotidis MI, Georgakilas AG. Molecular markers for cancer prognosis and treatment: have we struck gold? Cancer Lett. 2011;327(1):142–52.PubMedGoogle Scholar
- Modugno F, Ness RB, Chen C, Weiss NS. Inflammation and endometrial cancer: a hypothesis. Cancer Epidemiol Biomarkers Prev. 2005;14(12):2840–7.View ArticlePubMedGoogle Scholar
- Oertelt-Prigione S. The influence of sex and gender on the immune response. Autoimmun Rev. 2011, In Press, DOI:10.1016/j.autrev.2011.11.022.
- Remmers N, Bailey JM, Mohr AM, Hollingsworth MA. Molecular pathology of early pancreatic cancer. Cancer Biomark. 2011;9(1–6):421–40.Google Scholar
- Orsi NM, Gopichandran N, Ekbote UV, Walker JJ. Murine serum cytokines throughout the estrous cycle, pregnancy and post partum period. Anim Reprod Sci. 2006;96(1–2):54–65.View ArticlePubMedGoogle Scholar
- Rostene W, Dansereau MA, Godefroy D, Van Steenwinckel J, Reaux-Le Goazigo A, Melik-Parsadaniantz S, et al. Neurochemokines: a menage a trois providing new insights on the functions of chemokines in the central nervous system. J Neurochem. 2011;118(5):680–94.
- Baud V, Karin M. Signal transduction by tumor necrosis factor and its relatives. Trends Cell Biol. 2001;11(9):372–7.View ArticlePubMedGoogle Scholar
- Orsi NM, Tribe RM. Cytokine networks and the regulation of uterine function in pregnancy and parturition. J Neuroendocrinol. 2008;20(4):462–9.View ArticlePubMedGoogle Scholar
- Balkwill FR, Burke F. The cytokine network. Immunol Today. 1989;10(9):299–304.View ArticlePubMedGoogle Scholar
- Orsi NM. Cytokine networks in the establishment and maintenance of pregnancy. Hum Fertil (Camb). 2008;11(4):222–30.View ArticleGoogle Scholar
- Chaouat G, Ledee-Bataille N, Zourbas S, Ostojic S, Dubanchet S, Martal J, et al. Cytokines, implantation and early abortion: re-examining the Th1/Th2 paradigm leads to question the single pathway, single therapy concept. Am J Reprod Immunol. 2003;50(3):177–86.
- Barabasi AL, Gulbahce N, Loscalzo J. Network medicine: a network-based approach to human disease. Nat Rev Genet. 2011;12(1):56–68.PubMed CentralView ArticlePubMedGoogle Scholar
- Petricka JJ, Benfey PN. Reconstructing regulatory network transitions. Trends Cell Biol. 2011;21(8):442–51.PubMed CentralView ArticlePubMedGoogle Scholar
- Pe’er D. Bayesian network analysis of signaling networks: a primer. Sci STKE. 2005;2005(281):l4.Google Scholar
- Pe’er D, Hacohen N. Principles and strategies for developing network models in cancer. Cell. 2011;144(6):864–73.PubMed CentralView ArticlePubMedGoogle Scholar
- Grzegorczyk M. An introduction to Gaussian Bayesian networks. Methods Mol Biol. 2010;662:121–47.View ArticlePubMedGoogle Scholar
- Djebbari A, Quackenbush J. Seeded Bayesian Networks: constructing genetic networks from microarray data. BMC Syst Biol. 2008;2:57.PubMed CentralView ArticlePubMedGoogle Scholar
- Sachs K, Perez O, Pe’er D, Lauffenburger DA, Nolan GP. Causal protein-signaling networks derived from multiparameter single-cell data. Science. 2005;308(5721):523–9.View ArticlePubMedGoogle Scholar
- Ben-Gal I. Bayesian Networks. In: Ruggeri F, Faltin F, Kenett R, editors. Encyclopedia of Statistics in Quality & Reliability. 2013th ed. Chichester: Wiley and Sons; 2007.Google Scholar
- Heckerman D, Geiger D, Chickering D. Learning Bayesian networks: the combination of knowledge and statistical data. Mach Learn. 1995;20(95):197–234.Google Scholar
- Gyftodimos E, Flach P. Hierarchical Bayesian networks: a probabilistic reasoning model for structured domains. In: Proceedings of the ICML-2002 Workshop on Development of Representations 2002. 23–30.
- Friedman N, Linial M, Nachman I, Pe’er D. Using Bayesian networks to analyze expression data. J Comput Biol. 2000;7(3–4):601–20.View ArticlePubMedGoogle Scholar
- Husmeier D, Werhli AV. Bayesian integration of biological prior knowledge into the reconstruction of gene regulatory networks with Bayesian networks. Comput Syst Bioinformatics Conf. 2007;6:85–95.View ArticlePubMedGoogle Scholar
- Kuschner KW, Malyarenko DI, Cooke WE, Cazares LH, Semmes OJ, Tracy ER. A Bayesian network approach to feature selection in mass spectrometry data. BMC Bioinformatics. 2010;11:177.PubMed CentralView ArticlePubMedGoogle Scholar
- Aloy P, Russell RB. Structural systems biology: modelling protein interactions. Nat Rev Mol Cell Biol. 2006;7(3):188–97.View ArticlePubMedGoogle Scholar
- Buckley AR. Prolactin, a lymphocyte growth and survival factor. Lupus. 2001;10(10):684–90.View ArticlePubMedGoogle Scholar
- Reber PM. Prolactin and immunomodulation. Am J Med. 1993;95(6):637–44.View ArticlePubMedGoogle Scholar
- Cejkova P, Fojtikova M, Cerna M. Immunomodulatory role of prolactin in diabetes development. Autoimmun Rev. 2009;9(1):23–7.View ArticlePubMedGoogle Scholar
- Saha S, Gonzalez J, Rosenfeld G, Keiser H, Peeva E. Prolactin alters the mechanisms of B cell tolerance induction. Arthritis Rheum. 2009;60(6):1743–52.PubMed CentralView ArticlePubMedGoogle Scholar
- Losonsky GA, Ogra PL. Maternal-neonatal interactions and human breast milk. Prog Clin Biol Res. 1981;70:171–82.PubMedGoogle Scholar
- Brandtzaeg P. The secretory immunoglobulin system: regulation and biological significance. Focusing on human mammary glands. Adv Exp Med Biol. 2002;503:1–16.View ArticlePubMedGoogle Scholar
- Brandtzaeg P. Mucosal immunity: integration between mother and the breast-fed infant. Vaccine. 2003;21(24):3382–8.View ArticlePubMedGoogle Scholar
- Brandtzaeg P. The mucosal immune system and its integration with the mammary glands. J Pediatr. 2010;156(2 Suppl):S8–15.View ArticlePubMedGoogle Scholar
- Mestecky J. Homeostasis of the mucosal immune system: human milk and lactation. Adv Exp Med Biol. 2001;501:197–205.View ArticlePubMedGoogle Scholar
- Agarwal S, Karmaus W, Davis S, Gangur V. Immune markers in breast milk and fetal and maternal body fluids: a systematic review of perinatal concentrations. J Hum Lact. 2011;27(2):171–86.View ArticlePubMedGoogle Scholar
- Gao S, Wang X. Quantitative utilization of prior biological knowledge in the Bayesian network modeling of gene expression data. BMC Bioinformatics. 2011;12:359.PubMed CentralView ArticlePubMedGoogle Scholar
- Cagnin S, Biscuola M, Patuzzo C, Trabetti E, Pasquali A, Laveder P, et al. Reconstruction and functional analysis of altered molecular pathways in human atherosclerotic arteries. BMC Genomics. 2009;10:13.
- Hayes MP, Wang J, Norcross MA. Regulation of interleukin-12 expression in human monocytes: selective priming by interferon-gamma of lipopolysaccharide-inducible p35 and p40 genes. Blood. 1995;86(2):646–50.PubMedGoogle Scholar
- Roda G, Marocchi M, Sartini A, Roda E. Cytokine networks in ulcerative colitis. Ulcers. 2011, In Press.
- Arnaud L, Gorochov G, Charlotte F, Lvovschi V, Parizot C, Larsen M, et al. Systemic perturbation of cytokine and chemokine networks in Erdheim-Chester disease: a single-center series of 37 patients. Blood. 2011;117(10):2783–90.
- Bhavnani SK, Victor S, Calhoun WJ, Busse WW, Bleecker E, Castro M, et al. How cytokines co-occur across asthma patients: from bipartite network analysis to a molecular-based classification. J Biomed Inform. 2011;44 Suppl 1:S24–30.
- Horseman ND, Zhao W, Montecino-Rodriguez E, Tanaka M, Nakashima K, Engle SJ, et al. Defective mammopoiesis, but normal hematopoiesis, in mice with a targeted disruption of the prolactin gene. EMBO J. 1997;16(23):6926–35.
- Bouchard B, Ormandy CJ, Di Santo JP, Kelly PA. Immune system development and function in prolactin receptor-deficient mice. J Immunol. 1999;163(2):576–82.PubMedGoogle Scholar
- Penny LA. Monocyte chemoattractant protein 1 in luteolysis. Rev Reprod. 2000;5(2):63–6.View ArticlePubMedGoogle Scholar
- Wongdee K, Tulalamba W, Thongbunchoo J, Krishnamra N, Charoenphandhu N. Prolactin alters the mRNA expression of osteoblast-derived osteoclastogenic factors in osteoblast-like UMR106 cells. Mol Cell Biochem. 2011;349(1–2):195–204.View ArticlePubMedGoogle Scholar
- Armstrong DT, Miller LV, Knudsen KA. Regulation of lipid metabolism and progesterone production in rat corpora lutea and ovarian interstitial elements by prolactin and luteinizing hormone. Endocrinology. 1969;85(3):393–401.View ArticlePubMedGoogle Scholar
- Shipman LJ, Docherty AH, Knight CH, Wilde CJ. Metabolic adaptations in mouse mammary gland during a normal lactation cycle and in extended lactation. Q J Exp Physiol. 1987;72(3):303–11.View ArticlePubMedGoogle Scholar
- Wilde CJ, Knight CH, Flint DJ. Control of milk secretion and apoptosis during mammary involution. J Mammary Gland Biol Neoplasia. 1999;4(2):129–36.View ArticlePubMedGoogle Scholar
- Nishimori K, Young LJ, Guo Q, Wang Z, Insel TR, Matzuk MM. Oxytocin is required for nursing but is not essential for parturition or reproductive behavior. Proc Natl Acad Sci U S A. 1996;93(21):11699–704.PubMed CentralView ArticlePubMedGoogle Scholar
- Wagner KU, Young 3rd WS, Liu X, Ginns EI, Li M, Furth PA, et al. Oxytocin and milk removal are required for post-partum mammary-gland development. Genes Funct. 1997;1(4):233–44.
- Young 3rd WS, Shepard E, Amico J, Hennighausen L, Wagner KU, LaMarca ME, et al. Deficiency in mouse oxytocin prevents milk ejection, but not fertility or parturition. J Neuroendocrinol. 1996;8(11):847–53.
- Zicari A, Ticconi C, Realacci M, Cela O, Santangelo C, Pietropolli A, et al. Hormonal regulation of cytokine release by human fetal membranes at term gestation: effects of oxytocin, hydrocortisone and progesterone on tumour necrosis factor-alpha and transforming growth factor-beta 1 output. J Reprod Immunol. 2002;56(1–2):123–36.
- Maccio A, Madeddu C, Chessa P, Panzone F, Lissoni P, Mantovani G. Oxytocin both increases proliferative response of peripheral blood lymphomonocytes to phytohemagglutinin and reverses immunosuppressive estrogen activity. In Vivo. 2010;24(2):157–63.PubMedGoogle Scholar
- Amrani Y, Syed F, Huang C, Li K, Liu V, Jain D, et al. Expression and activation of the oxytocin receptor in airway smooth muscle cells: Regulation by TNFalpha and IL-13. Respir Res. 2010;11:104.
- Bluman EM, Bartynski KJ, Avalos BR, Caligiuri MA. Human natural killer cells produce abundant macrophage inflammatory protein-1 alpha in response to monocyte-derived cytokines. J Clin Invest. 1996;97(12):2722–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Hoh BL, Hosaka K, Downes DP, Nowicki KW, Fernandez CE, Batich CD, et al. Monocyte chemotactic protein-1 promotes inflammatory vascular repair of murine carotid aneurysms via a macrophage inflammatory protein-1alpha and macrophage inflammatory protein-2-dependent pathway. Circulation. 2011;124(20):2243–52.
- Goebeler M, Schnarr B, Toksoy A, Kunz M, Brocker EB, Duschl A, et al. Interleukin-13 selectively induces monocyte chemoattractant protein-1 synthesis and secretion by human endothelial cells. Involvement of IL-4R alpha and Stat6 phosphorylation. Immunology. 1997;91(3):450–7.
- Voest EE, Kenyon BM, O’Reilly MS, Truitt G, D’Amato RJ, Folkman J. Inhibition of angiogenesis in vivo by interleukin 12. J Natl Cancer Inst. 1995;87(8):581–6.View ArticlePubMedGoogle Scholar
- Liu YY, Slotine JJ, Barabasi AL. Controllability of complex networks. Nature. 2011;473(7346):167–73.View ArticlePubMedGoogle Scholar
- Yoshimoto T, Tsutsui H, Tominaga K, Hoshino K, Okamura H, Akira S, et al. IL-18, although antiallergic when administered with IL-12, stimulates IL-4 and histamine release by basophils. Proc Natl Acad Sci U S A. 1999;96(24):13962–6.
- Andrews AL, Nasir T, Bucchieri F, Holloway JW, Holgate ST, Davies DE. IL-13 receptor alpha 2: a regulator of IL-13 and IL-4 signal transduction in primary human fibroblasts. J Allergy Clin Immunol. 2006;118(4):858–65.View ArticlePubMedGoogle Scholar
- Arida A, Fragiadaki K, Giavri E, Sfikakis PP. Anti-TNF agents for Behcet’s disease: analysis of published data on 369 patients. Semin Arthritis Rheum. 2011;41(1):61–70.View ArticlePubMedGoogle Scholar
- Migliore A, Broccoli S, Bizzi E, Lagana B. Indirect comparison of the effects of anti-TNF biological agents in patients with ankylosing spondylitis by means of a mixed treatment comparison performed on efficacy data from published randomised, controlled trials. J Med Econ. 2012;15(3):473–80.View ArticlePubMedGoogle Scholar
- Benucci M, Saviola G, Baiardi P, Manfredi M, Sarzi Puttini P, Atzeni F. Determinants of risk infection during therapy with anti TNF-alpha blocking agents in rheumatoid arthritis. Open Rheumatol J. 2012;6:33–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Turkstra E, Ng SK, Scuffham PA. A mixed treatment comparison of the short-term efficacy of biologic disease modifying anti-rheumatic drugs in established rheumatoid arthritis. Curr Med Res Opin. 2011;27(10):1885–97.View ArticlePubMedGoogle Scholar
- Glabinski AR, Bielecki B, Kolodziejski P, Han Y, Selmaj K, Ransohoff RM. TNF-alpha microinjection upregulates chemokines and chemokine receptors in the central nervous system without inducing leukocyte infiltration. J Interferon Cytokine Res. 2003;23(8):457–66.View ArticlePubMedGoogle Scholar
- Meyer Zum Buschenfelde C, Cramer S, Trumpfheller C, Fleischer B, Frosch S. Trypanosoma cruzi induces strong IL-12 and IL-18 gene expression in vivo: correlation with interferon-gamma (IFN-gamma) production. Clin Exp Immunol. 1997;110(3):378–85.PubMed CentralView ArticlePubMedGoogle Scholar
- Rattenbacher B, Bohjanen PR. Evaluating posttranscriptional regulation of cytokine genes. Methods Mol Biol. 2012;820:71–89.PubMed CentralView ArticlePubMedGoogle Scholar
- Orsi NM, Ekbote UV, Walker JJ, Gopichandran N. Uterine and serum cytokine arrays in the mouse during estrus. Anim Reprod Sci. 2007;100(3–4):301–10.View ArticlePubMedGoogle Scholar
- Chandrashekar V, Bartke A, Wagner TE. Endogenous human growth hormone (GH) modulates the effect of gonadotropin-releasing hormone on pituitary function and the gonadotropin response to the negative feedback effect of testosterone in adult male transgenic mice bearing human GH gene. Endocrinology. 1988;123(6):2717–22.View ArticlePubMedGoogle Scholar
- Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Stat Soc B. 1995;57:289–300.Google Scholar
- Cooke EJ, Savage RS, Kirk PD, Darkins R, Wild DL. Bayesian hierarchical clustering for microarray time series data with replicates and outlier measurements. BMC Bioinformatics. 2011;12:399.PubMed CentralView ArticlePubMedGoogle Scholar
- Beal MJ, Falciani F, Ghahramani Z, Rangel C, Wild DL. A Bayesian approach to reconstructing genetic regulatory networks with hidden factors. Bioinformatics. 2005;21(3):349–56.View ArticlePubMedGoogle Scholar
- Jenssen TK, Laegreid A, Komorowski J, Hovig E. A literature network of human genes for high-throughput analysis of gene expression. Nat Genet. 2001;28(1):21–8.
- Hall M, Frank E, Holmes F, Pfahringer B, Reutemann P, Witten I. The WEKA data mining software: an update. SIGKDD Explorations. 2009;11(1):10–8.View ArticleGoogle Scholar
- Zhu J, Lum PY, Lamb J, GuhaThakurta D, Edwards SW, Thieringer R, et al. An integrative genomics approach to the reconstruction of gene networks in segregating populations. Cytogenet Genome Res. 2004;105(2–4):363–74.
- Friedman N, Goldszmidt M, Wyner A. Data analysis with Bayesian networks: a bootstrap approach. Proceedings of the Fifteenth conference on Uncertainty in artificial intelligence. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc; 1999:196-205. ISBN:1-55860-614-9.
- Russell SJ, Norvig P. Artificial Intelligence: A Modern Approach. 2nd ed. New Jersey: Prentice Hall; 2003.Google Scholar
- Shannon P, Markiel A, Ozier O, Baliga N, Wang J, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
- Bishop C. Pattern Recognition and Machine Learning. 2nd ed. Singapore: Springer Science and Business Media; 2007.Google Scholar
- Aldridge BB, Gaudet S, Lauffenburger DA, Sorger PK. Lyapunov exponents and phase diagrams reveal multi-factorial control over TRAIL-induced apoptosis. Mol Syst Biol. 2011;7:553.PubMed CentralView ArticlePubMedGoogle Scholar