 Research article
 Open Access
 Published:
Bayesian modeling suggests that IL12 (p40), IL13 and MCP1 drive murine cytokine networks in vivo
BMC Systems Biology volume 9, Article number: 76 (2015)
Abstract
Background
Cytokinehormone 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 cytokinehormone interactions in vivo using murine lactation as a dynamic, physiological model system.
Results
Circulatory levels of estrogen, progesterone, prolactin and twentythree 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), IL13 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 networkbased approaches in identifying novel cytokine pharmacological targets and in predicting the effects of their exogenous manipulation in inflammatory/immune disorders.
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 wellrecognized [1–10]. Many former reductionist studies have attempted to ascribe given physiological effects to the actions of one or small numbers of cytokines, but this strategy has met with conceptual difficulties stemming from the fact that their actions can be paradoxical at different concentrations and differ according to the prevailing hormonal milieu [3, 11]. Furthermore, the mechanistic insight offered by such studies in terms of relative cytokine interactions is also limited given that physiological responses are seldom governed by any one cytokine but rather by the combined influences of many. In this respect, these mediators are believed to operate as part of highly complex networks, wherein they exhibit antagonism, synergy and functional redundancy [12–14]. It is therefore clear that gaining a fuller understanding of cytokine function in any biological or clinical scenario rests with clarifying their interactions at a network level rather than relying on the increasingly inadequate T helper cell type 1/2 (Th1/Th2) paradigm [15].
Recent developments in highthroughput 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 [18–20]. 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.
A Bayesian network can be inferred from experimental data through the correlations between experimentallymeasured 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 priorknowledge graph ‘seed’ incorporating wellrecognized, literaturederived 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 [25–28].
While most Bayesianbased studies have focused on core biological processes using data from cell lines or tissuebased genomewide expression microarrays, less attention has been paid to physiological or clinically relevant in vivo systemicallyderived 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 hormonebased datasets: one drawn from lactating mice and the other from mice whose pups were removed at birth. With regard to the former network, murine lactation represents a unique model system to explore the dynamic physiological interactions between the hormonal environment and inflammation/immune function, with particular regard to the effects of prolactin (PRL), both a key driver of lactation and putative critical immunomodulator [30–33]. Since lactation PRL levels can conveniently be abrogated in vivo by removal of the suckling stimulus, this offers an excellent strategy to verify predictions made by in silico perturbation of the lactation network by comparing it with an independentlygenerated in vivo pupfree network with a physiologicallyinduced 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 perturbationdependent 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 [34–39]. The benefits of this approach included not having to use multiple pharmacological agents/antibodies whose effects would require titrating to concentration (likely in a nonlinear manner), and not being subject to nonphysiological 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.
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 IL1α, IL2, IL3, IL5, IL9, IL10, IL12 (p40), IL12 (p70), IL17, IFNγ, GCSF, GMCSF, KC, MCP1, MIP1α, MIP1β, RANTES, P_{4} and PRL (Table 1 and Additional file 1). Trends for IL13 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 IL6, IL9, IL12 (p70), IL13, GCSF, GMCSF, MCP1, MIP1α, MIP1β and TNFα. Furthermore, there was a marked trend for PRL, E_{2}, IL1α, IL2, IL9, IL10, IL12 (p40), IL12 (p70), IL13, eotaxin, GCSF, IFNγ, KC, MCP1, MIP1α and RANTES to decrease on day 2 of lactation. Most cytokine levels peaked at day 10 of lactation, particularly IL2 and MCP1, although the timing of this phenomenon differed for IL1α (days 16–24), IL9 (days 4–16), IL12 (p40) (day 24), KC (day 21) and RANTES (day 16). There followed a significant decrease (IL2, IFNγ, GCSF, GMCSF, MCP1, MIP1α, MIP1β; P < 0.05) or trend towards decreasing profiles (IL1β, IL3, IL4, IL5, IL6, IL10, IL12 (p40), IL12 (p70), IL13, IL17, eotaxin, KC, TNFα) on day 16. P_{4} 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 E_{2} and P_{4} were similar to those of NC animals.
Cytokine and hormone profiles during lactation cluster into three distinct timeseries
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 E_{2}, IL1β, IL9, IL12 (p40), P_{4} 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 (IL9, E_{2}, P_{4}, PRL) peaked around day 5, and tailed off steadily thereafter. Cluster 2 (IL1α, IL1β, IL12 (p40), IL12 (p70), IL17, MIP1β, RANTES) started off low and increased steadily to plateau from day 10 onwards, while Cluster 3 (eotaxin, GCSF, GMCSF, IL2, IL3, IL4, IL5, IL6, IL10, IL13, IFNγ, KC, MCP1, MIP1α, TNFα) behaved similarly but had a broad peak centered on day 10 of lactation.
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 IL3, E_{2} and eotaxin as the firstline parents and the other with IL12 (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γ, IL13, MCP1, MIP1α, MIP1β, 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 IL4 was orphaned. E_{2}, P_{4}, and PRL had a high probability of being present at elevated concentration relative to all other network components, with the exception of eotaxin and IL9 (Fig. 4).
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 overfitting 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 pseudodataset. 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 MIP1α) akin to those obtained from the seeded Bayesian models (Additional file 3). IL10, IL12 (p70) (a child of the seeded network parent IL12 (p40)), IL13, eotaxin and MIP1β were parents, and RANTES and TNFα were terminal nodes (with IL1α and IL2 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 Fscore 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.
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 P_{4} 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 IL13, MCP1 and IL2, 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 IL13 and MCP1), further perturbations were performed in order to explore the relative importance of other network elements.
IL13, MCP1 and IL12 (p40) are key driver nodes
As outlined above, IL13 and MCP1 were selected for perturbation by allocating them to a high concentration bin with a probability of 1 (based on their ~3 and 2fold increases in high concentration conditional probabilities, respectively, following in silico PRL depletion) (Table 3). They were perturbed both individually and in combination. Increasing IL13 concentration caused extensive network changes, including a shift in MCP1 and GMCSF 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 MIP1α, RANTES and TNFα became more medium. Note that the term ‘medium’ is used herein as shorthand to refer to the intermediate, midrange 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 MCP1 resulted in similarly extensive changes: among its children, GMCSF concentration increased markedly while MIP1β and RANTES became medium (Additional file 5 panel B). Further downstream effects included MIP1α concentration becoming high, KC moderately mediumhigh, and IFNγ and TNFα both more medium. PRL branch perturbation (combined PRL/IL13 and PRL/MCP1) 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/IL13/MCP1 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 IL3, E_{2} and eotaxin as firstline parents), IL12 (p40) was chosen for perturbation. IL12 (p40) perturbation in silico (by allocation to a high concentration bin) had a dramatic downstream impact: IL1β, IL12 (p70) and KC concentration increased, MIP1α became moderately high, and MIP1β, IFNγ and TNFα became medium (Additional file 7 panel A). Combined perturbations of both network branches (IL12 (p40) and eotaxin) also affected common downstream nodes (Additional file 8, panels A and B).
MCP1 can act synergistically with IL13 and/or IL12 (p40)
The combined perturbation of IL13 and MCP1 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 IL13 and MCP1, respectively). This suggests the presence of a synergistic interaction in which IL13 may coopt MCP1 given that it is also one of its parents (Table 3). However, more striking still were the synergistic effects of MCP1 and IL12 (p40) in relation to MIP1α (Fig. 5). Conditional probabilities in the MIP1α 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 MIP1β when IL12 (p40) and MCP1 were perturbed together, and for IFNγ when IL12 (p40) and IL13 were perturbed in combination (Additional file 9).
IL3 and PRL are potentially antagonistic
Similar perturbations of the parent node IL3 were also performed. Allocation of IL3 to low, medium and high concentration bins affected IL2, IL13 and MCP1 concentrations to some degree, although the most striking feature of these changes was that the greatest effects were noted when IL3 concentration was medium, rather than high or low (i.e. its effects were not linearly related to concentration; Additional file 10). Perturbing IL3 by allocating it to a high concentration bin resulted in an increase in IL13 concentration, as did perturbing PRL to a low concentration. This suggested that, as parents, IL3 and PRL may be antagonistic in terms of their effects on IL13, a notion supported by the intermediate conditional probabilities for IL13 concentration resulting from high levels of both these parents (Fig. 5).
MIP1β exhibits a biphasic response to eotaxin
Allocating eotaxin to a high concentration bin had minor effects, causing little more than a shift in MIP1β 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: IL9 concentration was reduced whereas MCP1 and MIP1β became more medium, an effect carried through downstream to GMCSF (Additional file 11 panel B). Intriguingly, the shift in MIP1β concentration was in the same direction, independent of whether eotaxin concentration was perturbed upwards or downwards, suggesting a concentrationdependent 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 lactationdependent rise in PRL). This physiological perturbation resulted in a fall in maternal serum concentrations of IL17 and a rise in KC on day 2 (corrected P < 0.05) (Fig. 6). By day 4 postpartum, the differences between females with and without pups were more pronounced: IL1α, IL12 (p40), IL17, IFNγ, GCSF, E_{2} and PRL levels were significantly higher in nursing dams (P < 0.05). Similar trends were also noted for IL2, IL5, IL9 and IL12 (p70).
Comparison of Bayesian networks
Lactation and pupfree networks share striking similarities
Strikingly, the pupfree dam Bayesian network retained the same core structural hubs as the lactation network (IFNγ, IL13, MCP1, MIP1α, MIP1β, and RANTES), with IL12 (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 (IL10), three orphan nodes (GCSF, IL4 and IL6) 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 pupfree dam networks were compared using accepted measures of network comparison (see Experimental Procedures), the Fscore 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.
As noted for the lactation network, the VBSSM variant for pupfree 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). IL12 (p40) and IL10 were principal parents, with IFNγ and MIP1β 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 Fcore 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 pupfree scenario, echoing the node status distribution in the in silico perturbed network driven by high concentration IL13 and MCP1 (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 (P_{4}, IL1α, IL2 and IL10 to low concentration bins). However, combined PRL/IL13/MCP1 perturbation correctly classified the qualitative (i.e. quantitatively high/medium/low concentrations relative to the equal frequency bins) nodal status of P_{4}, IL1α, IL10 (low concentration), RANTES, TNFα (medium), IFNγ (moderately high) and GMCSF (high). Close categorization was noted for MIP1β (predicted to be medium instead of moderately high), IL6 (moderately low instead of medium), IL17 (moderately low instead of low) and MIP1α (moderately high instead of high) (Table 3 and Fig. 7). In addition to having MCP1 as a parent, MIP1α, MIP1β and IFNγ were also children of IL12 (p40), and perturbing both of these parents by allocating them to a high concentration bin in silico led to their correct categorization in the pupfree network (Additional file 7 panel B and Fig. 7). Only the qualitative profile of GCSF, KC, IL5 and IL9 in the pupfree 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 pupfree 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 networkbased approach. This method revealed a robust lactation network structure featuring two main branches organized around principal parents (IL12 (p40), IL3, E_{2} and eotaxin), structural hubs (IFNγ, IL13, MCP1, MIP1α, MIP1β, 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], IL12 (p40) in those governing T cell and macrophage responses [42], IL13 in the pathophysiology of ulcerative colitis [43], IFNγ, IL12 (p40/p70) and MCP1 in ErdheimChester disease [44] and eotaxin, IFNγ, MIP1α and MIP1β in asthma [45]. A second network was generated from data drawn from pupfree 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 cytokinemediated 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 (IL12 (p40)) and two structural hubs (IL13 and MCP1) 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 P_{4} 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 P_{4}, only minor qualitative shifts in conditional probability values in downstream mediators (including those for IL13  PRL’s connection to the rest of the lactation network  and MCP1). 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 MCP1 are rather more unclear: in contrast to our findings, existing data implicate PRL as an inducer of MCP1 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 P_{4} supports PRL’s reported role in increasing serum P_{4} in rats [50].
In silico perturbation of downstream IL13, 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 IL13 input function was unable to achieve. It is tempting to speculate that oxytocin, a lactationrelated neurohypophyseal hormone which causes smooth muscle contraction during the letdown reflex, may represent an additional input given that the highest cytokine concentrations coincided with the period of maximal suckling/milk production [51–55]. This would be consistent with oxytocin’s known modulation of cytokine production and receptordependent activity in a range of settings [56, 57], and the ability of IL13 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 IL13 and MCP1 to high concentration bins in silico (akin to what occurred physiologically in the in vivo pupfree network) correctly (or closely) predicted the relative concentration status of 13 out of their 14 downstream cytokines in vivo. IL13 and MCP1 perturbation less accurately predicted MIP1α and MIP1β’s response to pup removal (based on their conditional probabilities). However, they were correctly categorized when their direct driver node parent (IL12 (p40) and MCP1) perturbations were used instead. These findings are consistent with the documented effects of MCP1 induction of MIP1α expression in murine aneurysm models and total IL12 induction of MIP1α in isolated human natural killer cells (albeit in the presence of IL15) [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 cellspecific 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 IL13 to be a selective inducer of MCP1 [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 IL13 and MCP1 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. IL12 (p40) and MCP1 perturbation independently increased the concentration of their direct child MIP1α, but this effect was markedly greater when these parents were perturbed together. Analogous, though less pronounced synergistic interactions were also noted for combined IL12 (p40) and IL13 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, IL3 and PRL had opposing actions on IL13, with IL3 seemingly increasing IL13 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 IL3), the effect on IL13 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, IL4 was peripheral to the main Bayesian network despite tracking IL13 levels, as established by the cluster analysis. This feature suggests functional redundancy, which is consistent with the fact that IL4 and IL13 are known to operate through the same receptor system [9, 14, 65]. It is worth noting that IL6 and GCSF were also orphaned in the pupfree 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 IL6 and GCSF 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 interrelationships, as demonstrated by the use of VBSSM models.
Another interesting phenomenon noted was the response of MIP1β to changes induced by in silico eotaxin manipulation. Intriguingly, perturbing eotaxin to either high or low concentrations both resulted in MIP1β concentration being higher. Concentrationdependent paradoxical effects are well recognized among cytokines and, in this instance, MIP1β responses point to a concentrationdependent 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 IL12 (p40), IL13 and MCP1 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 antibodybased 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 antiTNF therapy in the management of a range of autoimmune disorders [66–68]. Interestingly, comparative trials of various agents in rheumatoid arthritis suggest that agents such as these (e.g. etanercept, certolizumab) are superior to those targeting other cytokines such as IL1 (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 immunerelated side effects.
This Bayesian networkbased analysis has proved valuable in clarifying the complex structure and causal murine systemic cytokinehormone 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 contextspecific 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 cellspecific 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 MIP1β 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 IL13 and MCP1, between IL12 (p40) and IFNγ) [61, 71].
Conclusions
The identification of synergy, antagonism, functional redundancy and concentrationdependent 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 nonbiased 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 posttranscriptional 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) knockout 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 nonlactating mice (Experiment 2).
Animals
Experiment 1
Eight to tenweek 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 5565 % and 21.5 ± 1 °C. Females were naturally pairmated to 12–14 week old CD1 stud males of proven fertility following Whitten effectinduced 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 crossfostering 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 (timematched) (n = 8 in both groups; these 16 data points were used to construct pupfree networks). Seventysix 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 IL1α, IL1β, IL2, IL3, IL4, IL5, IL6, IL9, IL10, IL12 (p40), IL12 (p70), IL13, IL17, eotaxin, GCSF, GMCSF, IFNγ, KC, MCP1, MIP1α, MIP1β, RANTES and TNFα by multiplex immunoassay (BioRad Laboratories, Hemel Hempstead, Hertfordshire, UK) on a Luminex100 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 interplate 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 E_{2} and P_{4} were assayed by enzymelinked immunosorbent assay according to the manufacturer’s instructions (Alpha Diagnostic, San Antonio, Texas).
Data presentation and analysis
Data were expressed as pg/ml (cytokines, E_{2}) or ng/ml (PRL, P_{4}) ± SEM. All data distributions were assessed for normality by AndersonDarling tests. Basic analytical approaches were performed to highlight time courserelated changes in analyte profiles and to better appreciate the data before applying machine learning approaches. These were based on subsequent KruskallWallis/analysis of variance with post hoc Mann–WhitneyU/Fisher’s LSD tests, as appropriate. Pup removal data were similarly compared using ttests or Mann–WhitneyU 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, zscores 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 zscore P values (0.05 threshold) were corrected using the False Discovery Rate. Timeseries 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 timeseries were normalized to have zero mean and unit variance, the clustering analysis was sensitive only to their shape. A flexible, nonparametric 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, P_{1} and P_{2}, is associated with a conditional probability table containing probability of the state of P_{2} given the state of P_{1}. 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 P_{i} and Θ_{G} is the conditional distribution for each variable given its parents: P_{b}(P_{i}  P_{a}(P_{i})), where P_{a}(P_{i}) denotes the set of all parents of P_{i} 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, proteinprotein interaction databases, or any combination thereof, was used. The information incorporated in the seed network was solely restricted to wellestablished mousebased interactions, given that many cytokines exhibit speciesspecific 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/productsservices/pharmalifesciences/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 structurerelated factor in the scoring metric is given in Heckerman et al. (1995) [24]. Prior to performing the Bayesian network analysis, zscore normalizations were applied to the raw data in Matlab. A machine learning algorithm (implemented in the WEKAbased opensource 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 highconfidence 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 zscore 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 pupfree) following the independent categorization described above for each dataset.
Timedependent 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 binbased concentration states than the firstorder 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 oneoff 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 edgespecific 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 nonroot 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 nonparametric bootstrapping was applied (100 operations) to address potential overfitting in the Bayesian analysis, wherein multiple data sets were created by resampling 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 beliefnetwork 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 statespace 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 statespace models to help reverseengineer interactions between proteins or transcriptional networks from time series data, thereby explicitly modelling the progression of a system over a multistep time course. In statespace 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 oneshot 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 firstorder 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 zscore). In this regard, the zscore 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 1p. Accordingly, a larger zscore 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 zscore 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, IL3, IL12 (p40), IL13, MCP1 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, IL13 and MCP1 were perturbed (given a high concentration status), both individually and in combination. As a final step, IL3 and IL12 (p40) were perturbed alone and in combination with IL13, MCP1 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 pupremoval (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.rproject.org/web/packages/catnet/index.html), a categorical Bayesian network inference framework was used to systematically assess the structural similarities between the lactation and pupfree networks by computing the Fscore, which represents the harmonic average of specificity and sensitivity, and accounts for interstructural 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:
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 c_{1}, c_{2}, .., c_{k} has complexity c_{1} x c_{2} x … x c_{k}. 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 pupfree 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).
Abbreviations
 BHC:

Bayesian hierarchical clustering
 E_{2} :

estrogen
 GCSF:

granulocytecolony stimulating factor
 GMCSF:

granulocyte macrophagecolony stimulating factor
 IL:

interleukin
 IFN:

interferon
 KC:

keratinocyte chemoattractant
 MCP:

monocyte chemoattractant protein
 MIP:

macrophage inflammatory protein
 P_{4} :

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
References
 1.
Crivellato E. The role of angiogenic growth factors in organogenesis. Int J Dev Biol. 2011;55(4–5):365–75.
 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.
LopezPedrera C, PerezSanchez C, RamosCasals M, SantosGonzalez M, RodriguezAriza A, Cuadrado MJ. Cardiovascular risk in systemic autoimmune diseases: epigenetic mechanisms of immune regulatory functions. Clin Dev Immunol. 2012;2012:9746–8.
 4.
Ma S, Ma CC. Recent development in pleiotropic effects of statins on cardiovascular disease through regulation of transforming growth factorbeta superfamily. Cytokine Growth Factor Rev. 2011;22(3):167–75.
 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.
 6.
Modugno F, Ness RB, Chen C, Weiss NS. Inflammation and endometrial cancer: a hypothesis. Cancer Epidemiol Biomarkers Prev. 2005;14(12):2840–7.
 7.
OerteltPrigione 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.
 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.
 10.
Rostene W, Dansereau MA, Godefroy D, Van Steenwinckel J, ReauxLe Goazigo A, MelikParsadaniantz 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.
 12.
Orsi NM, Tribe RM. Cytokine networks and the regulation of uterine function in pregnancy and parturition. J Neuroendocrinol. 2008;20(4):462–9.
 13.
Balkwill FR, Burke F. The cytokine network. Immunol Today. 1989;10(9):299–304.
 14.
Orsi NM. Cytokine networks in the establishment and maintenance of pregnancy. Hum Fertil (Camb). 2008;11(4):222–30.
 15.
Chaouat G, LedeeBataille N, Zourbas S, Ostojic S, Dubanchet S, Martal J, et al. Cytokines, implantation and early abortion: reexamining 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 networkbased approach to human disease. Nat Rev Genet. 2011;12(1):56–68.
 17.
Petricka JJ, Benfey PN. Reconstructing regulatory network transitions. Trends Cell Biol. 2011;21(8):442–51.
 18.
Pe’er D. Bayesian network analysis of signaling networks: a primer. Sci STKE. 2005;2005(281):l4.
 19.
Pe’er D, Hacohen N. Principles and strategies for developing network models in cancer. Cell. 2011;144(6):864–73.
 20.
Grzegorczyk M. An introduction to Gaussian Bayesian networks. Methods Mol Biol. 2010;662:121–47.
 21.
Djebbari A, Quackenbush J. Seeded Bayesian Networks: constructing genetic networks from microarray data. BMC Syst Biol. 2008;2:57.
 22.
Sachs K, Perez O, Pe’er D, Lauffenburger DA, Nolan GP. Causal proteinsignaling networks derived from multiparameter singlecell data. Science. 2005;308(5721):523–9.
 23.
BenGal I. Bayesian Networks. In: Ruggeri F, Faltin F, Kenett R, editors. Encyclopedia of Statistics in Quality & Reliability. 2013th ed. Chichester: Wiley and Sons; 2007.
 24.
Heckerman D, Geiger D, Chickering D. Learning Bayesian networks: the combination of knowledge and statistical data. Mach Learn. 1995;20(95):197–234.
 25.
Gyftodimos E, Flach P. Hierarchical Bayesian networks: a probabilistic reasoning model for structured domains. In: Proceedings of the ICML2002 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.
 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.
 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.
 29.
Aloy P, Russell RB. Structural systems biology: modelling protein interactions. Nat Rev Mol Cell Biol. 2006;7(3):188–97.
 30.
Buckley AR. Prolactin, a lymphocyte growth and survival factor. Lupus. 2001;10(10):684–90.
 31.
Reber PM. Prolactin and immunomodulation. Am J Med. 1993;95(6):637–44.
 32.
Cejkova P, Fojtikova M, Cerna M. Immunomodulatory role of prolactin in diabetes development. Autoimmun Rev. 2009;9(1):23–7.
 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.
 34.
Losonsky GA, Ogra PL. Maternalneonatal interactions and human breast milk. Prog Clin Biol Res. 1981;70:171–82.
 35.
Brandtzaeg P. The secretory immunoglobulin system: regulation and biological significance. Focusing on human mammary glands. Adv Exp Med Biol. 2002;503:1–16.
 36.
Brandtzaeg P. Mucosal immunity: integration between mother and the breastfed infant. Vaccine. 2003;21(24):3382–8.
 37.
Brandtzaeg P. The mucosal immune system and its integration with the mammary glands. J Pediatr. 2010;156(2 Suppl):S8–15.
 38.
Mestecky J. Homeostasis of the mucosal immune system: human milk and lactation. Adv Exp Med Biol. 2001;501:197–205.
 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.
 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.
 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 interleukin12 expression in human monocytes: selective priming by interferongamma of lipopolysaccharideinducible p35 and p40 genes. Blood. 1995;86(2):646–50.
 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 ErdheimChester disease: a singlecenter 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 cooccur across asthma patients: from bipartite network analysis to a molecularbased classification. J Biomed Inform. 2011;44 Suppl 1:S24–30.
 46.
Horseman ND, Zhao W, MontecinoRodriguez 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 receptordeficient mice. J Immunol. 1999;163(2):576–82.
 48.
Penny LA. Monocyte chemoattractant protein 1 in luteolysis. Rev Reprod. 2000;5(2):63–6.
 49.
Wongdee K, Tulalamba W, Thongbunchoo J, Krishnamra N, Charoenphandhu N. Prolactin alters the mRNA expression of osteoblastderived osteoclastogenic factors in osteoblastlike UMR106 cells. Mol Cell Biochem. 2011;349(1–2):195–204.
 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.
 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.
 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.
 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.
 54.
Wagner KU, Young 3rd WS, Liu X, Ginns EI, Li M, Furth PA, et al. Oxytocin and milk removal are required for postpartum mammarygland 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 factoralpha and transforming growth factorbeta 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.
 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 IL13. Respir Res. 2010;11:104.
 59.
Bluman EM, Bartynski KJ, Avalos BR, Caligiuri MA. Human natural killer cells produce abundant macrophage inflammatory protein1 alpha in response to monocytederived cytokines. J Clin Invest. 1996;97(12):2722–7.
 60.
Hoh BL, Hosaka K, Downes DP, Nowicki KW, Fernandez CE, Batich CD, et al. Monocyte chemotactic protein1 promotes inflammatory vascular repair of murine carotid aneurysms via a macrophage inflammatory protein1alpha and macrophage inflammatory protein2dependent pathway. Circulation. 2011;124(20):2243–52.
 61.
Goebeler M, Schnarr B, Toksoy A, Kunz M, Brocker EB, Duschl A, et al. Interleukin13 selectively induces monocyte chemoattractant protein1 synthesis and secretion by human endothelial cells. Involvement of IL4R 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.
 63.
Liu YY, Slotine JJ, Barabasi AL. Controllability of complex networks. Nature. 2011;473(7346):167–73.
 64.
Yoshimoto T, Tsutsui H, Tominaga K, Hoshino K, Okamura H, Akira S, et al. IL18, although antiallergic when administered with IL12, stimulates IL4 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. IL13 receptor alpha 2: a regulator of IL13 and IL4 signal transduction in primary human fibroblasts. J Allergy Clin Immunol. 2006;118(4):858–65.
 66.
Arida A, Fragiadaki K, Giavri E, Sfikakis PP. AntiTNF agents for Behcet’s disease: analysis of published data on 369 patients. Semin Arthritis Rheum. 2011;41(1):61–70.
 67.
Migliore A, Broccoli S, Bizzi E, Lagana B. Indirect comparison of the effects of antiTNF 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.
 68.
Benucci M, Saviola G, Baiardi P, Manfredi M, Sarzi Puttini P, Atzeni F. Determinants of risk infection during therapy with anti TNFalpha blocking agents in rheumatoid arthritis. Open Rheumatol J. 2012;6:33–7.
 69.
Turkstra E, Ng SK, Scuffham PA. A mixed treatment comparison of the shortterm efficacy of biologic disease modifying antirheumatic drugs in established rheumatoid arthritis. Curr Med Res Opin. 2011;27(10):1885–97.
 70.
Glabinski AR, Bielecki B, Kolodziejski P, Han Y, Selmaj K, Ransohoff RM. TNFalpha microinjection upregulates chemokines and chemokine receptors in the central nervous system without inducing leukocyte infiltration. J Interferon Cytokine Res. 2003;23(8):457–66.
 71.
Meyer Zum Buschenfelde C, Cramer S, Trumpfheller C, Fleischer B, Frosch S. Trypanosoma cruzi induces strong IL12 and IL18 gene expression in vivo: correlation with interferongamma (IFNgamma) production. Clin Exp Immunol. 1997;110(3):378–85.
 72.
Rattenbacher B, Bohjanen PR. Evaluating posttranscriptional regulation of cytokine genes. Methods Mol Biol. 2012;820:71–89.
 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.
 74.
Chandrashekar V, Bartke A, Wagner TE. Endogenous human growth hormone (GH) modulates the effect of gonadotropinreleasing 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.
 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.
 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.
 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.
 78.
Jenssen TK, Laegreid A, Komorowski J, Hovig E. A literature network of human genes for highthroughput 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.
 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:196205. ISBN:1558606149.
 82.
Russell SJ, Norvig P. Artificial Intelligence: A Modern Approach. 2nd ed. New Jersey: Prentice Hall; 2003.
 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.
 85.
Aldridge BB, Gaudet S, Lauffenburger DA, Sorger PK. Lyapunov exponents and phase diagrams reveal multifactorial control over TRAILinduced apoptosis. Mol Syst Biol. 2011;7:553.
Acknowledgements
Invaluable advice was provided by Dr Ian Overton, Dr Rebecca Ward, Dr Megha Padi, Dr Benjamin HaibeKains, 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.
Author information
Additional information
Competing interests
The author(s) declare that they have no competing interests.
Authors’ contributions
NMO and JG designed and supervised the project. HM, MC, NMO and SLF conducted animal work, experimental analyses and/or basic statistical analyses. RSS performed the cluster analysis while TD and JA constructed the Bayesian networks and related statistical analyses. MC, NMO, SLF and TD interpreted the Bayesian networks. JG, MC, NMO, RSS, SLF and TD contributed to writing the manuscript. All authors read and approved the final manuscript.
Sarah L. Field and Tathagata Dasgupta contributed equally to this work.
Additional files
Additional file 1:
Changes in cytokine and hormone profiles throughout lactation; Cytokine, E _{ 2 } (pg/ml), PRL and P _{ 4 } (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. (TIF 644 kb)
Additional file 2:
Prior network utilised to develop the Bayesian network. The prior network was learned from the literature relating to murine cytokine interactions only. Nodes are displayed within circles, with edges representing known interactions. Due to the nature of the discovery of the prior network, edges do not represent directionality. (TIF 1064 kb)
Additional file 3:
Circular variational Bayesian statespace model (VBSSM) of cytokine interactions during physiological lactation; Although it was created entirely from the experimental data with no inferences drawn from prior knowledge, much of the core nodal structure present in the seeded network also emerges in this acyclic graph (A). Circular VBSSM of cytokine interactions in the pupfree group. The nodes p40 and p70 refer to the IL12 (p40) subunit and (p70) heterodimer, respectively (B). (TIF 914 kb)
Additional file 4:
Bayesian lactation network perturbation by deterministically decreasing PRL concentration; Includes conditional probability indicators for each component node (color coding is as described for Fig. 3 in the main text). Bars beside each node represent conditional probabilities (low to high, from right to left). (TIF 3620 kb)
Additional file 5:
Bayesian lactation network perturbation by deterministically increasing IL13 concentration (A), MCP1 concentration (B) and in combination (C). (TIF 3966 kb)
Additional file 6:
Bayesian lactation network perturbation by PRL branch (combined PRL/IL13/MCP1) perturbation. (TIF 6170 kb)
Additional file 7:
Bayesian lactation network perturbation by deterministically increasing IL12 (p40) concentration alone (A) or in conjunction with increasing MCP1 (B). (TIF 2996 kb)
Additional file 8:
Bayesian lactation network perturbation by deterministically increasing IL12 (p40) concentration in the presence of decreased (A) or increased (B) eotaxin. (TIF 3060 kb)
Additional file 9:
Changes in conditional probability associated with perturbation of eotaxin and IL12 (p40) as parent nodes relative to those induced by perturbing MCP1 and IL13 as hubs. Subscripts indicate the direction of the perturbation (−low, +high). (DOCX 14 kb)
Additional file 10:
Conditional probabilities associated with IL3 perturbation. Subscripts indicate the direction of the perturbation (Llow, MMedium, Hhigh). (DOCX 14 kb)
Additional file 11:
Bayesian lactation network perturbation by deterministically increasing (A) and decreasing (B) eotaxin concentration. (DOCX 14 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Field, S.L., Dasgupta, T., Cummings, M. et al. Bayesian modeling suggests that IL12 (p40), IL13 and MCP1 drive murine cytokine networks in vivo . BMC Syst Biol 9, 76 (2015). https://doi.org/10.1186/s1291801502263
Received:
Accepted:
Published:
Keywords
 Bayesian network
 Cytokine
 Hormone
 Lactation
 Mouse