Open Access

Bayesian modeling suggests that IL-12 (p40), IL-13 and MCP-1 drive murine cytokine networks in vivo

  • Sarah L. Field1,
  • Tathagata Dasgupta2,
  • Michele Cummings1,
  • Richard S. Savage3,
  • Julius Adebayo2, 4,
  • Hema McSara1,
  • Jeremy Gunawardena2 and
  • Nicolas M. Orsi1Email author
Contributed equally
BMC Systems Biology20159:76

DOI: 10.1186/s12918-015-0226-3

Received: 29 October 2015

Accepted: 31 October 2015

Published: 9 November 2015

Abstract

Background

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.

Results

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).

Conclusion

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.

Keywords

Bayesian network Cytokine Hormone Lactation Mouse

Background

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 [110]. 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 [1214]. 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 [15].

Recent developments in high-throughput analytical platforms have revolutionized the quality and quantity of data available from in vivo experiments. The large number of analytes measurable in single samples provides the opportunity to explore their interrelationships in both physiological and pathophysiological processes [16, 17]. In this regard, Bayesian networks provide an attractive methodology for analyzing such complex biological data [1820]. Given that many biologists are unlikely to be familiar with probabilistic graphical models, a word of introduction to Bayesian networks is warranted. A Bayesian network is a directed acyclic graph whose nodes are the variables of interest (herein, a cytokine/hormone), each of which can have a range of quantitative values which are typically discretized into a small number of bins, such as ‘low’, ‘medium’ and ‘high’. The directed edges in the graph (represented as arrows between nodes) reflect likely causal relationships between nodes. The nature of these causal relationships is captured by the graph’s underlying conditional probability table (CPT) which details the probabilities for any given node to fall into each of the different (in this case, concentration) bins given the status of its parent nodes (i.e. those directly upstream). The underlying CPT does not change upon perturbation; rather the marginal probability of that node displaying a certain behavior changes (Fig. 1). The illustrative Bayesian network shown in Fig. 1 has five variables (vertices/“nodes”). Node E is not causally influenced by any of the others, nor does it causally influence them, so this node has no edges entering or leaving it (i.e. it is ‘orphaned’). By contrast, nodes B and C are solely influenced by A, so each has a single connecting edge. However, they respond to A in quite different ways. Based on the conditional probability tables, if the value of A is categorized into a low concentration bin, then B has a marginal probability of 0.8 of falling into a high concentration bin, while C has a probability of 0.75 of falling into a low one. The status of D is influenced by both B and C and accordingly has two incoming edges. This approach offers the scope for an intricate description of D’s behavior, based on the conditional probabilities associated with the allocation of its data to high or low concentration bins which, in turn, depends on the state of A, B and C. The conditional probability tables for each node represent relative (rather than absolute) concentrations.
Fig. 1

Illustrative Bayesian network describing causal relationships between five variables, with their associated conditional probability tables. The values of each variable have been discretized into low (l), medium (m) and high (h) bins. The notation P(B = l|A) refers to the probability of B being in the low value bin, conditional upon the value of A (which itself can be l, m or h). Note that with one parent only (i.e. the case for nodes B and C), both row and column probabilities sum to 1, whereas with multiple inputs (i.e. in the case of D), only the rows sum to 1. The nodes are colored green (high), white (medium) and red (low) to illustrate in silico perturbation where A falls into a high concentration bin (probability 0.9) and E falls into a low concentration bin (probability 0.8). The implications of this are demonstrated through the conditional probability tables associated with each downstream node wherein having D in a given concentration state is dictated by the particular combination of states of its parents B and C as dictated by their corresponding overall (marginal) probabilities (the entries as captured in the conditional probability tables do not change on intervention), as summarized in the histogram attached to each node. Note that the overall marginal probability of D being in a particular state (which changes under intervention) is the sum as summarized in the histogram attached to each node. Readers more familiar with Bayesian networks will note that the sub-network structures (B → A, A → C, B → D, C → D) and (C → A, A → B, B → D, C → D) are both Markov equivalent with the present one (A → B, A → C, B → D, C → D). When we condition on A, the other nodes B and C become independent no matter whether A is a tail-to-head or tail-to-tail intermediate node (this scenario differs from that seen with head-to-head node D). In other words, these three sub-networks specify the same independent assumptions belonging to the same equivalence class and the true causal network can possibly be any one of the sub-network solutions. However, causal networks based on observation alone (i.e. without intervention, which is an important tool for inferring causality) can be still partially constructed. In the present study, prior knowledge seed network edge directionalities were assigned by a modified depth-first search algorithm which helped to choose the sub-network from an equivalence class containing more than one Markov equivalent member as suggested here (see later)

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) [21]. 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 [24]. 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 [2528].

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 [3033]. 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).

Results

We chose to characterize the physiological putative changes in, and interactions between, systemic hormone and cytokine concentrations using murine lactation as a model, since this system offered the benefit of featuring both intimate and highly variable physiological immunoendocrine interactions [3439]. The benefits of this approach included not having to use multiple pharmacological agents/antibodies whose effects would require titrating to concentration (likely in a non-linear manner), and not being subject to non-physiological interactions or inducing global changes unrelated to the system of interest. Similarly, our approach deliberately aimed to avoid using multiple knockout models on both pragmatic and functional grounds given the compensatory redundancy of cytokine networks and/or the unknown postpartum traits of such models. Our design was based on two independent data sets: one covering cytokine and hormonal profiles over 7 time points throughout lactation in order to generate a Bayesian model of mediator interactions, and another abrogated lactation data set which provided a biological validation platform for assessing the predictive power of the former network in determining the physiological profile changes elicited by pup removal. The methodological workflow is described in Fig. 2; please see the methods section for experimental details.
Fig. 2

Experimental methodology utilized in the present study; initial biological analysis (via fluid-phase multiplex immunoassay, ELISA and radioimmunoassay) was followed by development and exploration of Bayesian networks. In vivo physiological perturbation was performed to validate the in silico networks.

Biological analysis

Most circulatory cytokines and hormones vary markedly throughout lactation

Significant changes (P < 0.05, following correction for multiple comparisons) in concentration during lactation were noted for IL-1α, IL-2, IL-3, IL-5, IL-9, IL-10, IL-12 (p40), IL-12 (p70), IL-17, IFN-γ, G-CSF, GM-CSF, KC, MCP-1, MIP-1α, MIP-1β, RANTES, P4 and PRL (Table 1 and Additional file 1). Trends for IL-13 were also noted (P = 0.075). Cytokine levels showed a broad tendency to be increased on day 1 of lactation relative to naturally cycling (NC) concentrations, and significant increases were observed in IL-6, IL-9, IL-12 (p70), IL-13, G-CSF, GM-CSF, MCP-1, MIP-1α, MIP-1β and TNF-α. Furthermore, there was a marked trend for PRL, E2, IL-1α, IL-2, IL-9, IL-10, IL-12 (p40), IL-12 (p70), IL-13, eotaxin, G-CSF, IFN-γ, KC, MCP-1, MIP-1α and RANTES to decrease on day 2 of lactation. Most cytokine levels peaked at day 10 of lactation, particularly IL-2 and MCP-1, although the timing of this phenomenon differed for IL-1α (days 16–24), IL-9 (days 4–16), IL-12 (p40) (day 24), KC (day 21) and RANTES (day 16). There followed a significant decrease (IL-2, IFN-γ, G-CSF, GM-CSF, MCP-1, MIP-1α, MIP-1β; P < 0.05) or trend towards decreasing profiles (IL-1β, IL-3, IL-4, IL-5, IL-6, IL-10, IL-12 (p40), IL-12 (p70), IL-13, IL-17, eotaxin, KC, TNF-α) on day 16. P4 concentrations rose in line with those of PRL, peaked on day 4, fell to very low levels by day 16, and increased towards weaning. Weaning concentrations for E2 and P4 were similar to those of NC animals.
Table 1

Cytokine and hormone concentrations throughout lactation

Mediator

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)

IL-1α

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

IL-1β

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

IL-2

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

IL-3

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

IL-4

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

IL-5

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

IL-6

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

IL-9

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

IL-10

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

IL-12 (p40)

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

IL-12 (p70)

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

IL-13

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

IL-17

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

Eotaxin

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

G-CSF

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

GM-CSF

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

IFN-γ

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

KC

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

MCP-1

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

MIP-1α

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

MIP-1β

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

RANTES

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

TNF-α

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

E2

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

P4

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

PRL

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, E2 (pg/ml), PRL and P4 (ng/ml) concentrations throughout lactation (mean ± SEM). Groups that do not share a common superscript letter are significantly different from each other (P < 0.05); NC - naturally cycling

Cytokine and hormone profiles during lactation cluster into three distinct time-series

In order to reveal the temporal structure in the data, we examined correlations between analytes and clustered the time series. Significant correlations were noted across the array of mediators investigated, except for E2, IL-1β, IL-9, IL-12 (p40), P4 and PRL (Fig. 3, panel a). These relationships were used to inform the cluster analysis (Fig. 3, panel b) which revealed that analyte profiles fell into three clusters: Cluster 1 (IL-9, E2, P4, PRL) peaked around day 5, and tailed off steadily thereafter. Cluster 2 (IL-1α, IL-1β, IL-12 (p40), IL-12 (p70), IL-17, MIP-1β, RANTES) started off low and increased steadily to plateau from day 10 onwards, while Cluster 3 (eotaxin, G-CSF, GM-CSF, IL-2, IL-3, IL-4, IL-5, IL-6, IL-10, IL-13, IFN-γ, KC, MCP-1, MIP-1α, TNF-α) behaved similarly but had a broad peak centered on day 10 of lactation.
Fig. 3

Correlations between cytokine and hormone mediators; Panel a Heat map demonstrating the correlations between cytokines, E2, P4 and PRL. Red coloring indicates positive correlations, while blue coloring indicates negative correlations. The strength of color reflects the strength of correlation (as indicated by the scale) – numbers in brackets refer to cluster identities. Panel b clusters indicating relationships between mediators. Three clusters were identified: Cluster 1 - IL-9, E2, P4, PRL; Cluster 2 - IL-1α, IL-1β, IL-12 (p40), IL-12 (p70), IL-17, MIP-1β, RANTES; and cluster 3 - eotaxin, G-CSF, GM-CSF, IL-2, IL-3, IL-4, IL-5, IL-6, IL-10, IL-13, IFN-γ, KC, MCP-1, MIP-1α, TNF-α

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

A Bayesian network from the lactation dataset was generated incorporating a prior knowledge network (Additional file 2) with experimental data, wherein node (i.e. mediator) interactions were portrayed as directed edges implying likely causal relationships. Each node was associated with a set of conditional probabilities which determined its status (i.e. probable relative concentration) dependent upon the status of its parents. This network organized into two main branches, one with IL-3, E2 and eotaxin as the first-line parents and the other with IL-12 (p40) as the principal parent node (Fig. 4). The network itself was assembled around six structural hubs (i.e. possible signal integrators; defined as nodes with >1 input and output edges totaling ≥5) comprising IFN-γ, IL-13, MCP-1, MIP-1α, MIP-1β, and RANTES. The terminal node was TNF-α, which was connected - both directly and indirectly - to each of these hubs. A total of 42 directed edges (35 of which were of high confidence) representing cytokine causal relationships connected all but one node: only IL-4 was orphaned. E2, P4, and PRL had a high probability of being present at elevated concentration relative to all other network components, with the exception of eotaxin and IL-9 (Fig. 4).
Fig. 4

Bayesian network showing cytokine interrelationships in murine lactation. The nodes are color-coded according to the conditional probability of corresponding mediator relative concentrations being high (green), low (red) or medium (white) concentration given the state(s) of their parent nodes; however, within the present network no nodes fell into the medium category. Relative to the white color, the normalized concentration (low or high) determines the intensity of the node color. Very high confidence level edges (causal connecting lines between nodes) are colored in green (or grey if lower than the 0.9 cutoff), based upon the confidence analysis of the Bayesian result

There were more samples (53) than nodes (26, initially) in the present study. This compares well with previous biological studies. Djebbari and Quackenbush (2008) [21] learned a Bayesian network with 63 features from 38 samples, while Gao and Wang (2011) [40] developed a Bayesian network learned with 36 features from 46 samples. Ideally, the use of many more samples would have allowed us to create the most robust model, but this need had to be offset against making an ethical and appropriate use of animals. Nonetheless, the highly stringent bootstrapping approach used herein ensured robustness of the model. Its purpose was to resample the data with replacement in order to avoid over-fitting so as to address the issue of using a limited number of samples. The final network contained 25 nodes wherein the maximum theoretical number of parents for each node could be 24. Networks were learned for each generated pseudo-dataset. Model averaging of these was then performed to obtain edge confidence scores. The overall bootstrap confidence was estimated by evaluating how many times relative to the total number of iterations a particular feature of interest (i.e. directed edge, Markov relation, etc.) appeared. The final network had an overall confidence of 0.9 (i.e. the learned network features were present in at least 90 networks out of 100 iterations). Increasing the confidence threshold from 0.7 to 0.9 had no effect, confirming the result with high confidence (i.e. strongly supported by the data). This was further confirmed by the fact that a Variational Bayesian State Space Model (VBSSM) network (which did not account for prior knowledge of known mediator interactions in their construction) revealed conservation of the core network structure, including hubs (except MIP-1α) akin to those obtained from the seeded Bayesian models (Additional file 3). IL-10, IL-12 (p70) (a child of the seeded network parent IL-12 (p40)), IL-13, eotaxin and MIP-1β were parents, and RANTES and TNF-α were terminal nodes (with IL-1α and IL-2 as peripheral termini). These topological changes were expected given that the VBSSM network only featured a subset of the original network nodes following automated construction (Additional file 3). The inference process was performed without prior knowledge bias of mediator interactions. As indicated in Table 2, the F-score value, high number of true negatives, and the relatively high sensitivity value when compared to the seeded network in Fig. 4 show that there is a relatively high topographical similarity between the two. This is further confirmed by the retention of major regulatory hubs and the fact that all the edges associated with TNF-α are incoming.
Table 2

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

Specificity

0.94

0.95

0.88

Sensitivity

0.30

0.14

0.17

F-Score

0.46

0.25

0.28

True positives

13

6

5

False positives

18

14

15

False negatives

30

36

25

True negatives

263

268

108

Bayesian networks model conditional independence so that accurately removing arcs from all possible connections is an important measurement for accuracy. This is reflected in the specificity, which is close to one if true negatives (TNs) are high and false positives (FPs) are low. For the purpose of internal validation (in addition to the experimental aspects), the VBSSM-based results obtained without prior knowledge were compared to those from seeded Bayesian learning. It is important to note that VBSSM results were derived under a strict confidence level check. However, we cannot expect a high agreement in true positives (TPs) for the structural comparisons performed for both the lactation and pup-free data but the high TN and specificity values are encouraging (first two data columns). The last column represents the network structural comparison within the VBSSM analysis between lactation and pup-free data

Exploration of networks

Perturbation of PRL in silico

PRL perturbation in silico in the lactation model was chosen as a starting point given that this hormone is widely considered as having potent immunomodulatory properties. Perturbation was achieved by allocating its low concentration bin a conditional probability of 1. This perturbation revealed that, apart from a fall in P4 profile, the downstream nodes, hubs and edges remained qualitatively unchanged (i.e. no/minimal directional shift (up or down) in conditional probabilities such that nodal status also remained unchanged (Table 3), an observation supported by the lack of correlation between PRL and the majority of other mediators (Fig. 3). The terms ‘increased’ and ‘decreased’ are used as less cumbersome terminology throughout this manuscript to describe changes in node conditional probabilities following perturbation (e.g. ‘decreased’ effectively relates to a shift in node conditional probability towards a low concentration bin). The only exceptions were IL-13, MCP-1 and IL-2, whose concentrations ‘increased’, indicating a degree of negative regulation of these cytokines by PRL (Table 3 and Additional file 4). Due to the relatively minor influence of in silico PRL depletion on the overall network (except, most notably, its negative regulation of the immediate downstream hubs IL-13 and MCP-1), further perturbations were performed in order to explore the relative importance of other network elements.
Table 3

Changes in conditional probability associated with perturbation of PRL and structural hubs

Cytokine

Low

Medium

High

PRL

0.065

0.692

0.243

P4

0.066

0.786

0.148

P4-PRL

0.778

0.111

0.111

P4-PRL+IL-13

0.778

0.111

0.111

P4-PRL+MCP-1

0.778

0.111

0.111

P4-PRL+IL-13+MCP-1

0.778

0.111

0.111

IL-13

0.400

0.499

0.101

IL-13-PRL

0.379

0.348

0.272

IL-13-PRL+IL-13

0.000

0.000

1.000

IL-13-PRL+MCP-1

0.379

0.348

0.272

IL-13-PRL+IL-13+MCP-1

0.000

0.000

1.000

IL-2

0.274

0.602

0.124

IL-2-PRL

0.274

0.481

0.246

IL-2-PRL+IL-13

0.111

0.111

0.778

IL-2-PRL+MCP-1

0.274

0.481

0.246

IL-2-PRL+IL-13+MCP-1

0.111

0.111

0.778

MCP-1

0.431

0.458

0.111

MCP-1-PRL

0.418

0.386

0.196

MCP-1-PRL+IL-13

0.200

0.229

0.571

MCP-1-PRL+MCP-1

0.000

0.000

1.000

MCP-1-PRL+IL-13+MCP-1

0.000

0.000

1.000

GM-CSF

0.359

0.526

0.115

GM-CSF-PRL

0.353

0.469

0.178

GM-CSF-PRL+IL-13

0.225

0.316

0.458

GM-CSF-PRL+MCP-1

0.111

0.111

0.778

GM-CSF-PRL+IL-13+MCP-1

0.111

0.111

0.778

MIP-1α

0.235

0.687

0.078

MIP-1α-PRL

0.250

0.644

0.106

MIP-1α-PRL+IL-13

0.250

0.521

0.229

MIP-1α-PRL+MCP-1

0.270

0.360

0.370

MIP-1α-PRL+IL-13+MCP-1

0.270

0.360

0.370

MIP-1β

0.607

0.208

0.186

MIP-1β-PRL

0.570

0.224

0.205

MIP-1β-PRL+IL-13

0.434

0.293

0.273

MIP-1β-PRL+MCP-1

0.321

0.357

0.321

MIP-1β-PRL+IL-13+MCP-1

0.321

0.357

0.321

KC

0.165

0.706

0.129

KC-PRL

0.168

0.692

0.140

KC-PRL+IL-13

0.183

0.635

0.181

KC-PRL+MCP-1

0.201

0.569

0.230

KC-PRL+IL-13+MCP-1

0.201

0.569

0.230

IL-1α

0.626

0.223

0.151

IL-1α-PRL

0.605

0.230

0.165

IL-1α-PRL+IL-13

0.528

0.256

0.216

IL-1α-PRL+MCP-1

0.465

0.281

0.254

IL-1α-PRL+IL-13+MCP-1

0.465

0.281

0.254

IL-10

0.648

0.230

0.122

IL-10-PRL

0.622

0.245

0.133

IL-10-PRL+IL-13

0.526

0.301

0.173

IL-10-PRL+MCP-1

0.448

0.348

0.204

IL-10-PRL+IL-13+MCP-1

0.448

0.348

0.204

IFN-γ

0.504

0.299

0.197

IFN-γ-PRL

0.465

0.299

0.236

IFN-γ-PRL+IL-13

0.325

0.337

0.338

IFN-γ-PRL+MCP-1

0.352

0.325

0.323

IFN-γ-PRL+IL-13+MCP-1

0.316

0.341

0.343

IL-5

0.517

0.327

0.156

IL-5-PRL

0.487

0.336

0.177

IL-5-PRL+IL-13

0.377

0.386

0.237

IL-5-PRL+MCP-1

0.398

0.375

0.228

IL-5-PRL+IL-13+MCP-1

0.370

0.390

0.240

IL-17

0.513

0.391

0.096

IL-17-PRL

0.486

0.406

0.108

IL-17-PRL+IL-13

0.390

0.467

0.144

IL-17-PRL+MCP-1

0.407

0.454

0.138

IL-17-PRL+IL-13+MCP-1

0.384

0.470

0.146

IL-6

0.608

0.306

0.086

IL-6-PRL

0.585

0.319

0.096

IL-6-PRL+IL-13

0.506

0.371

0.123

IL-6-PRL+MCP-1

0.521

0.361

0.119

IL-6-PRL+IL-13+MCP-1

0.501

0.374

0.124

RANTES

0.483

0.276

0.241

RANTES-PRL

0.454

0.287

0.259

RANTES-PRL+IL-13

0.373

0.322

0.305

RANTES-PRL+MCP-1

0.328

0.347

0.326

RANTES-PRL+IL-13+MCP-1

0.327

0.347

0.326

TNF-α

0.372

0.317

0.311

TNF-α-PRL

0.364

0.321

0.316

TNF-α-PRL+IL-13

0.346

0.328

0.325

TNF-α-PRL+MCP-1

0.341

0.331

0.329

TNF-α-PRL+IL-13+MCP-1

0.340

0.331

0.329

Conditional probabilities are given to each of three bins (low, medium and high). Highest conditional probability values in any given bin indicate a greater likelihood of the relevant mediator’s concentration being in that concentration bin. The large font denotes the node of interest; small font indicates the perturbed parent node; a ‘-’ sign indicates perturbation by reducing concentration while a ‘+’ sign indicates perturbation by increasing concentration

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)

The combined perturbation of IL-13 and MCP-1 on a background of low PRL had particularly marked effects on IFN-γ, resulting in an increase in its high concentration bin conditional probability (0.343), an effect greater than that achieved by each parent perturbation in isolation (0.338 and 0.323 from 0.236 for IL-13 and MCP-1, respectively). This suggests the presence of a synergistic interaction in which IL-13 may co-opt MCP-1 given that it is also one of its parents (Table 3). However, more striking still were the synergistic effects of MCP-1 and IL-12 (p40) in relation to MIP-1α (Fig. 5). Conditional probabilities in the MIP-1α high concentration bin were much greater when both were combined (0.363 and 0.370 compared to 0.600 in combination). Comparable, though less striking effects were noted for MIP-1β when IL-12 (p40) and MCP-1 were perturbed together, and for IFN-γ when IL-12 (p40) and IL-13 were perturbed in combination (Additional file 9).
Fig. 5

Synergistic and antagonistic relationships between mediators. Schematic representation of a MCP-1 and IL-12 (p40) synergy on MIP-1α with associated conditional probability table, and b antagonism between PRL and IL-3 on IL-13

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

In order to determine experimentally the impact of lactation on cytokine and hormone profiles relative to a baseline, samples were obtained from dams whose pups were removed at birth (i.e. from animals which would not have exhibited a lactation-dependent rise in PRL). This physiological perturbation resulted in a fall in maternal serum concentrations of IL-17 and a rise in KC on day 2 (corrected P < 0.05) (Fig. 6). By day 4 post-partum, the differences between females with and without pups were more pronounced: IL-1α, IL-12 (p40), IL-17, IFN-γ, G-CSF, E2 and PRL levels were significantly higher in nursing dams (P < 0.05). Similar trends were also noted for IL-2, IL-5, IL-9 and IL-12 (p70).
Fig. 6

Systemic cytokine profiles following parturition with and without lactation. Comparison of cytokine and PRL/P4/E2 concentrations in mice with/without pups on a day 2 and b day 4 of lactation (*P < 0.05)

Comparison of Bayesian networks

Lactation and pup-free networks share striking similarities

Strikingly, the pup-free dam Bayesian network retained the same core structural hubs as the lactation network (IFN-γ, IL-13, MCP-1, MIP-1α, MIP-1β, and RANTES), with IL-12 (p40) as the principal parent and TNF-α as the terminal node, despite expected differences in network topology (Fig. 7). This second network featured an additional parent (IL-10), three orphan nodes (G-CSF, IL-4 and IL-6) unconnected to the main network and a total of 42 edges (32 of high confidence) which connected 23 of the 26 nodes. When the lactation and pup-free dam networks were compared using accepted measures of network comparison (see Experimental Procedures), the F-score was 0.861 and the total complexities of these networks were 379 and 375, respectively. These scores imply a striking similarity in terms of topology and complexity between both networks despite their being generated from entirely independent data sets.
Fig. 7

Bayesian network displaying cytokine interrelationships associated with pup-free dams in vivo. Color coding of both nodes and edges is as described for Fig. 4. Symbols indicate if the nodal concentration bin was correctly categorized by in silico perturbation. - Correctly/closely categorized by PRL perturbation alone. - Correctly/closely categorized by PRL perturbation in combination with IL-12 (p40), IL-13 and MCP-1, or by any other single mediator perturbation. X – Not correctly categorized by any perturbation. * - Could not be evaluated due to being upstream of the perturbations performed

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.

Discussion

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 [41], IL-12 (p40) in those governing T cell and macrophage responses [42], IL-13 in the pathophysiology of ulcerative colitis [43], IFN-γ, IL-12 (p40/p70) and MCP-1 in Erdheim-Chester disease [44] and eotaxin, IFN-γ, MIP-1α and MIP-1β in asthma [45]. 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 [30]. 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 [50].

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 [5155]. 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 [58].

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 [61] 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 [64]), 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 [6668]. 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) [69], 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 [25]. 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 [70]. 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].

Conclusions

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) [72]. 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.

Methods

Methodology

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).

Animals

Experiment 1

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).

Experiment 2

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.

Ethics statement

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 [73], 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 [9]. 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 [74], 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 [75]. 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 [76]. 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 [18]. 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 [77]. 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) [24]. 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 [78]; 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 [81]. Several metrics and search algorithms were employed in combination in order to optimize the sensitivity of the results to the learning procedure scoring metric [82]; 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 [24]. 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) [83].

VBSSM model generation

In order to further validate the robustness of this model, a VBSSM was built in Matlab [77]. 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 [84], 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 [77]. 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 [63]. 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

In silico network perturbations preserved network structure as they only propagated the disturbance along a fixed network topology. By contrast, in vivo perturbation (i.e. pup removal) resulted in network topological changes, as anticipated given that the second network was constructed from an independent data set. The R package Catnet (http://cran.r-project.org/web/packages/catnet/index.html), a categorical Bayesian network inference framework was used to systematically assess the structural similarities between the lactation and pup-free networks by computing the F-score, which represents the harmonic average of specificity and sensitivity, and accounts for inter-structural edge appearances/disappearances. It is expressed in terms of the number of true positive (TP), true negative (TN), false positive (FP) and false negative (FN) edges, and represents a statistical score of similarity between any two networks in terms of edge connection features. In this instance, TP, TN, FP and FN refer to the directed edges of the in vivo perturbed network with respect to the lactation network (taken as the standard of truth). The following formula was used:
$$ \mathrm{F}\hbox{-} \mathrm{score} = 2\times \mathrm{specificity}\times \mathrm{sensitivity}/\left(\mathrm{specificity} + \mathrm{sensitivity}\right) $$

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 [85], 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).

Notes

Abbreviations

BHC: 

Bayesian hierarchical clustering

E2

estrogen

G-CSF: 

granulocyte-colony stimulating factor

GM-CSF: 

granulocyte macrophage-colony stimulating factor

IL: 

interleukin

IFN: 

interferon

KC: 

keratinocyte chemoattractant

MCP: 

monocyte chemoattractant protein

MIP: 

macrophage inflammatory protein

P4

progesterone

PRL: 

prolactin

RANTES: 

regulated upon activation, normal T cell expressed and secreted

SE: 

standard error

SEM: 

standard error of the mean

TNF: 

tumor necrosis factor

VBSSM: 

Variational Bayesian state space model

Declarations

Acknowledgements

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.

Authors’ Affiliations

(1)
Women’s Health Research Group, Section of Cancer & Pathology, Leeds Institute of Cancer & Pathology, St James’s University Hospital
(2)
Department of Systems Biology, Harvard Medical School
(3)
Systems Biology Centre, University of Warwick
(4)
IDeA Labs, Department of Computer Science, TMCB 1124, Brigham Young University

References

  1. Crivellato E. The role of angiogenic growth factors in organogenesis. Int J Dev Biol. 2011;55(4–5):365–75.View ArticlePubMedGoogle Scholar
  2. 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.
  3. 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
  4. 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
  5. 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
  6. 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
  7. 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.
  8. Remmers N, Bailey JM, Mohr AM, Hollingsworth MA. Molecular pathology of early pancreatic cancer. Cancer Biomark. 2011;9(1–6):421–40.Google Scholar
  9. 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
  10. 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.
  11. Baud V, Karin M. Signal transduction by tumor necrosis factor and its relatives. Trends Cell Biol. 2001;11(9):372–7.View ArticlePubMedGoogle Scholar
  12. 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
  13. Balkwill FR, Burke F. The cytokine network. Immunol Today. 1989;10(9):299–304.View ArticlePubMedGoogle Scholar
  14. Orsi NM. Cytokine networks in the establishment and maintenance of pregnancy. Hum Fertil (Camb). 2008;11(4):222–30.View ArticleGoogle Scholar
  15. 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.
  16. 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
  17. Petricka JJ, Benfey PN. Reconstructing regulatory network transitions. Trends Cell Biol. 2011;21(8):442–51.PubMed CentralView ArticlePubMedGoogle Scholar
  18. Pe’er D. Bayesian network analysis of signaling networks: a primer. Sci STKE. 2005;2005(281):l4.Google Scholar
  19. Pe’er D, Hacohen N. Principles and strategies for developing network models in cancer. Cell. 2011;144(6):864–73.PubMed CentralView ArticlePubMedGoogle Scholar
  20. Grzegorczyk M. An introduction to Gaussian Bayesian networks. Methods Mol Biol. 2010;662:121–47.View ArticlePubMedGoogle Scholar
  21. Djebbari A, Quackenbush J. Seeded Bayesian Networks: constructing genetic networks from microarray data. BMC Syst Biol. 2008;2:57.PubMed CentralView ArticlePubMedGoogle Scholar
  22. 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
  23. 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
  24. 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
  25. 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.
  26. 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
  27. 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
  28. 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
  29. Aloy P, Russell RB. Structural systems biology: modelling protein interactions. Nat Rev Mol Cell Biol. 2006;7(3):188–97.View ArticlePubMedGoogle Scholar
  30. Buckley AR. Prolactin, a lymphocyte growth and survival factor. Lupus. 2001;10(10):684–90.View ArticlePubMedGoogle Scholar
  31. Reber PM. Prolactin and immunomodulation. Am J Med. 1993;95(6):637–44.View ArticlePubMedGoogle Scholar
  32. Cejkova P, Fojtikova M, Cerna M. Immunomodulatory role of prolactin in diabetes development. Autoimmun Rev. 2009;9(1):23–7.View ArticlePubMedGoogle Scholar
  33. 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
  34. Losonsky GA, Ogra PL. Maternal-neonatal interactions and human breast milk. Prog Clin Biol Res. 1981;70:171–82.PubMedGoogle Scholar
  35. 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
  36. Brandtzaeg P. Mucosal immunity: integration between mother and the breast-fed infant. Vaccine. 2003;21(24):3382–8.View ArticlePubMedGoogle Scholar
  37. Brandtzaeg P. The mucosal immune system and its integration with the mammary glands. J Pediatr. 2010;156(2 Suppl):S8–15.View ArticlePubMedGoogle Scholar
  38. Mestecky J. Homeostasis of the mucosal immune system: human milk and lactation. Adv Exp Med Biol. 2001;501:197–205.View ArticlePubMedGoogle Scholar
  39. 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
  40. 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
  41. 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.
  42. 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
  43. Roda G, Marocchi M, Sartini A, Roda E. Cytokine networks in ulcerative colitis. Ulcers. 2011, In Press.
  44. 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.
  45. 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.
  46. 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.
  47. 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
  48. Penny LA. Monocyte chemoattractant protein 1 in luteolysis. Rev Reprod. 2000;5(2):63–6.View ArticlePubMedGoogle Scholar
  49. 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
  50. 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
  51. 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
  52. 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
  53. 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
  54. 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.
  55. 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.
  56. 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.
  57. 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
  58. 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.
  59. 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
  60. 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.
  61. 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.
  62. 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
  63. Liu YY, Slotine JJ, Barabasi AL. Controllability of complex networks. Nature. 2011;473(7346):167–73.View ArticlePubMedGoogle Scholar
  64. 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.
  65. 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
  66. 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
  67. 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
  68. 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
  69. 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
  70. 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
  71. 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
  72. Rattenbacher B, Bohjanen PR. Evaluating posttranscriptional regulation of cytokine genes. Methods Mol Biol. 2012;820:71–89.PubMed CentralView ArticlePubMedGoogle Scholar
  73. 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
  74. 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
  75. 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
  76. 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
  77. 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
  78. 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.
  79. 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
  80. 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.
  81. 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.
  82. Russell SJ, Norvig P. Artificial Intelligence: A Modern Approach. 2nd ed. New Jersey: Prentice Hall; 2003.Google Scholar
  83. 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.
  84. Bishop C. Pattern Recognition and Machine Learning. 2nd ed. Singapore: Springer Science and Business Media; 2007.Google Scholar
  85. 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

Copyright

© Field et al. 2015