|Home | About | Journals | Submit | Contact Us | Français|
Pathogenic bacteria such as Bordetella bronchiseptica modulate host immune responses to enable their establishment and persistence; however, the immune response is generally successful in clearing these bacteria. Here, we model the dynamic outcome of the interplay between host immune components and B. bronchiseptica virulence factors. The model extends our previously published interaction network of B. bronchiseptica and includes the existing experimental information on the relative timing of IL10 and IFNγ activation in the form of qualitative inequalities. The current model improves the previous one in two directions: (i) by augmenting the network with new nodes with specific function in T helper cell differentiation and effector mechanisms and (ii) by using a dynamic approach that allows us to quantify node states and mechanisms revealed to be important from our previous model. The model makes predictions about the time scales of each process, the activity thresholds of each node and novel regulatory interactions. For example, the model predicts that the activity threshold of IL4 is higher than that of IL12 and that pro-inflammatory cytokines regulate the activity of Th2 cells. Some of these predictions are supported by the literature, and many can serve as targets of future experiments.
The regulation of immune responses is a complex system of mechanisms, which has rarely been holistically explored. Immunological studies reveal the abundance of selected immune components at a few time points but relatively less effort is being made to understand and quantify the regulation among these components. Network modelling can assist in this process by integrating the behaviour of multiple components and addressing questions that are not yet accessible to experimental analysis. In our previous work (Thakar et al. 2007), we assembled networks of interactions between a mammalian immune system and two Gram-negative respiratory pathogens (Bordetella bronchiseptica and B. pertussis) and formulated discrete dynamic models that reproduce the basic features of the infection such as clearance or persistence of bacteria. Here, we extend and further parametrize the B. bronchiseptica model to achieve a more quantitative agreement with the available experimental data. Our study is motivated by the current state of immunology in which qualitative and comparative data are more abundant than quantitative data such as rate constants and by the increasing evidence that qualitative information can provide mechanistic knowledge and constructive hypotheses. We employ a hybrid dynamic modelling approach that incorporates combinatorial regulation as well as continuous degradation of immune components and bacterial effector functions. We impose constraints to select outcomes of the dynamic model based on known observations of bacterial and cytokine time courses and show that the model reproduces complex responses surprisingly well.
Using an infectious agent to study the regulation of immune responses not only defines the signal that initiates the response but also provides a variety of constraints that can be applied to the dynamic model. Bacteria persist within their hosts by subverting phagocytosis by immune cells, interfering with antigen processing or presentation, or by promoting anti-inflammatory or immunosuppressive responses that normally function to terminate the protective effector immune responses of the host (Mills 2004). We used B. bronchiseptica because there is a fair amount of information available about the immune response to this pathogen. Bordetella bronchiseptica colonizes the respiratory tracts of its hosts, adhering to ciliated epithelia and spreading via respiratory droplets. Bordetella bronchiseptica naturally infects wild and domesticated animals including mice (Cameron et al. 1980; van den Akker 1998) and can cause persistent disease typified by atrophic rhinitis in pigs and kennel cough in dogs. Certain aspects of this infection, such as the effect of the type III secretion system (TTSS) of B. bronchiseptica, have been well characterized. We used the time-course data in Pilione & Harvill (2006) to design constraints for the analysis of regulatory elements.
A network describing the immune response to B. bronchiseptica during its infection of the lower respiratory system was assembled in Thakar et al. (2007) using the literature of the mammalian immune system (Janeway et al. 2001) and specific immune responses to B. bronchiseptica. The overall topology of the network is maintained here. We have improved (figure 1) the spatial description of the components to address local and systemic responses (by the identification of two compartments; see text 1 and figure 1 in the electronic supplementary material) and the T helper cell regulatory network (by the incorporation of two additional experimentally measured nodes and several new edges; see text 2 in the electronic supplementary material). Compartment I represents the site of the infection (the lungs) and compartment II corresponds to the sites of T0 cell and B cell activation (see text 1 in the electronic supplementary material). Bacteria and the components of the immune system (i.e. immune cells and cytokines) are represented as network nodes, and interactions, regulatory relationships between components and their transport across compartments are represented as directed edges starting from the source (regulator) node and ending in the target node. We incorporated regulatory relationships that modulate a process (or an unspecified mediator of a process) as edges directed towards another edge. The regulatory effect of each edge is classified into activation or inhibition, and is represented by an incoming arrow (black) or an incoming blunt segment (red), respectively.
We assumed that the bacteria activating the chain of immune responses shown in figure 1 express the generic adhesins and toxins of the sequenced Bordetella species, and did not assign independent nodes for these virulence factors. Among species-specific virulence factors, we included O-antigen and TTSS as separate nodes because their specific functions have been characterized. Bordetella bronchiseptica can enter lymph nodes and lymphoid tissues (Gueirard et al. 2003), and TTSS is functional both in and outside of the lungs (Lee & Schneewind 1999), hence we indicate the presence of TTSS in both compartments.
The resulting 32-node network is described here, illustrating the sequence of processes activated after bacterial invasion (figure 1). The first immune mechanisms that respond to B. bronchiseptica's lipopolyssacharide (LPS) are receptor-mediated activation of respiratory epithelial cells and the alternative complement pathway. Receptor-mediated signalling activates the production of cytokines and chemokines. Many of these, including IL1, IL6, TNFα and TNFβ, are pro-inflammatory and recruit polymorphonuclear leucocytes (PMNs) to the site of infection (Mann et al. 2004). Activation of PMNs produces cytokines, which, in turn, recruit more phagocytes.
Following the commencement of phagocytosis, dendritic cells (the antigen-presenting cells in the network) are activated by cytokines and are directed to the nearest lymph node (compartment II) where they present antigens to T0 cells (naive T cells). Signalling pathways activated in both T0 and dendritic cells during this interaction induce the production of cytokines. The context-dependent cytokine profile directs the differentiation of T0 cells into either T helper type 1 (Th1) cells or T helper type 2 (Th2) cells. In this study, we use the following representative cytokines: IL12 and IL4 required during differentiation of T0 cells in compartment II (lymph nodes), and IFNγ and IL10 produced at the site of infection (lungs) by Th1 and Th2 cells, respectively. We analysed a variety of possible interactions among these cytokines based on the literature (refer to text 2 in the electronic supplementary material). This analysis predicts mutual inhibition between IL12 and IL4, and suggests two ways of inhibition of IFNγ by IL10, first, in synergy with TTSS as previously suggested (Pilione & Harvill 2006) and, second, by an interaction that is co-regulated by the other Th2-related cytokines (represented by IL4). Note that the dichotomy of differentiation into a Th1 or Th2 cell is at the single-cell level and the population-level effects are often described as the ratio of Th1 and Th2 cells. Differentiated T helper cells are then transported to the lungs where they engage in effector mechanisms.
The TTSS has been hypothesized to favour the generation of T regulatory cells (Tr; McGuirk et al. 2002; Skinner et al. 2005). During bacterial infections, Tr cells regulate Th1 responses especially via the production of IL10. Additionally, IL10 is an anti-inflammatory cytokine, and inhibits the production of pro-inflammatory cytokines. Th2 cells activate B cells in compartment II, which then proliferate and produce antibodies against bacterial antigens. Opsonization of bacteria by complement-fixing antibodies leads to the activation of the classical complement pathway. PMNs and macrophages (MP) express receptors that recognize complement-coated bacteria and bacterial antigens bound to antibodies; both interactions activate phagocytes. Pro-inflammatory cytokines and IFNγ attract more phagocytes, namely macrophages and natural killer (NK) cells, to the site of infection, where they are activated and eliminate antibody-bound bacteria. Thus, the Th2- and Th1-mediated antigen-specific responses constitute a positive feedback to the recruitment and activation of phagocytes, resulting in the elimination of the pathogen. The main mechanism of natural clearance of the bordetellae is via phagocytosis by activated phagocytes (Pishko et al. 2004; Kirimanjeswara et al. 2005).
The B. bronchiseptica LPS contains long repeats of O-antigen, which inhibit the activation of the alternative complement pathway (Burns et al. 2003). The B. bronchiseptica TTSS plays several roles. First, it induces the necrosis of PMNs (Yuk et al. 2000; Stockbauer et al. 2003). Second, the TTSS inhibits IFNγ production via the activation of IL10 during the first two weeks of the infection (Pilione & Harvill 2006), thus delaying the differentiation of T0 cells into Th1 cells and facilitating their differentiation into Th2 cells (Skinner et al. 2005). Third, TTSS has been shown to modulate dendritic cell function (Skinner et al. 2005). We used published experimental datasets (Pilione & Harvill 2006), in which the cytokine production in infections in mice with wild-type (WT) and TTSS-defective strains of B. bronchiseptica is studied.
We integrate the static network (figure 1) with time-course data for concentrations of IL10, IFNγ and bacterial numbers to develop a time-dependent dynamic model. The experiments consist of inoculation of mice with 5×105 colony forming units (CFU) of either WT B. bronchiseptica or a mutant derivative that lacks the TTSS due to an engineered deletion of the bscN gene encoding the ATPase required for the TTSS. In the rest of this paper, ‘mutant’ refers to the TTSS-defective strain unless otherwise specified. After inoculation, the lungs were excised to determine bacterial numbers and spleens were used to determine cytokine concentrations induced in response to B. bronchiseptica. The mutant bacteria are defective in persistence compared with the WT strain, suggesting that the inhibition of IFNγ by TTSS may be part of a strategy to survive inside the host.
We used the datasets of bacterial growth (reproduced as figure 2a) and the corresponding IL10 and IFNγ concentrations (reproduced as figure 2 in the electronic supplementary material) in the case of infection by WT and mutant bacteria. We could not aim for a complete quantitative agreement between experiments and model because the experimental data points of the time courses are not regularly dispersed, the cytokine and bacterial time courses originate from different tissues and the experimentally measured components represent the collective response of many signalling pathways. Thus, the model parameters were constrained to reproduce the qualitative features of the experimental bacterial growth curves and cytokine time courses.
Two constraints, bacterial clearance and B, Th1 and Th2 cell activation at least once during the simulations, were required to be satisfied for both WT and mutant simulations based on general knowledge of the B. bronchiseptica infection time course.
The simulated cytokine time courses in compartment I were required to satisfy the following conditions. First, there should be an association between bacterial clearance and high IFNγ concentration. Second, in the WT infection, the IL10 activation time point should be earlier than the IFNγ activation time point; the opposite order of activation should be observed in the mutant simulation. Third, the mutant simulation should reproduce the peculiar qualitative behaviour of IFNγ in which the concentration of IFNγ increases, followed by a drop and a second increase in concentration. All simulations satisfying these conditions were also visually screened for the IFNγ, IL10 and bacterial time courses.
The following qualitative features were assessed for the simulated bacterial growth curves in compartment I (figure 2a): in the WT infection, first, an initial exponential increase in bacterial numbers; second, a plateau or slow decrease following the initial increase indicating slow phagocytosis; and third, exponential decrease in the bacterial numbers following the plateau. In the mutant infection, first, an exponential increase in bacterial numbers, second, the absence of a plateau or a very small plateau compared with the WT bacterial growth curve and, third, faster decrease in the bacterial concentrations after the peak, leading to earlier clearance of mutant bacteria compared with the WT bacteria. For convenience, and following the classification into three infection phases introduced in Thakar et al. (2007), we refer to the behaviour leading to the first, second and third conditions as phases1 I, II and III, respectively. We selected parameter sets in which WT bacteria were cleared 1.8–2.8 times slower than mutant bacteria. The range was used because the experimental observations are based on a very small sample size (three mice) and are done sparsely. The lower limit of 1.8 (50/28) was based on the observation that WT bacteria clearance may be as early as day 50, while the 2.8 (70/24.5) upper limit was based on the estimate that mutant bacteria are cleared between days 21 and 28 (averaged to 24.5).
The conditions regarding the qualitative features of the bacterial growth curves were used to decide the functional form of terms in the differential equations describing the bacterial numbers and the strength of phagocytosis, and to parametrize these equations. All other conditions were used to filter out spurious parameter combinations.
During the course of infection, immune components and bacteria/virulence factors have a time-dependent expression. Given the scarcity of kinetic and quantitative characterization of the processes involved in the interaction network (figure 1), we employed a hybrid method that provides a bridge between discrete and continuous approaches. This method is especially suited to the problem because the model in Thakar et al. (2007) indicates that the decay and/or desensitization of immune components is crucial for their dynamics. The hybrid method allows us to describe node concentrations and activities by using coarse-grained parameters instead of molecular-level kinetic or binding parameters, which are unknown. As we will see, although the parameters are at a higher level than the kinetics of elementary reactions, they describe the biological mechanisms in sufficient detail to reveal their interdependency. Each node is described with two variables, a two-valued (Boolean) discrete variable Xi that stands for the activity of the node and a continuous variable that stands for the concentration of the node. Following Glass (1975), the time evolution of the continuous variables is described by piecewise linear differential equations that combine logical (Boolean) rules for the activation of the nodes with linear (free) decay (Glass 1975; De Jong et al. 2004; Chaves et al. 2006). The general form of the equation is , where Bi is a Boolean function that depends on (a subset of) the discrete variables in the network and γi is the decay rate of . The Boolean functions and decay rates are individually defined for each node. The differential equations of all nodes are given in table 1. At each instant t, the discrete variable Xi is defined by the continuous variable according to the threshold rule
In the case of cytokines, components expected to be produced faster than others, we augmented the general equation by another scaling factor α, first described in Chaves et al. (2006), that multiplies both terms in the differential equation. The parameters θi, γi and αi were analysed to find out their effective range and how they reflect the regulation of the immune response.
We described the continuous variables corresponding to bacteria (denoted ) and phagocytosis (denoted ) by ordinary differential equations (ODEs). The bacterial concentration is a dimensionless variable in the unit range and follows . The discrete variable defined for bacteria expresses whether the bacterial concentration is above the threshold level θBT to activate epithelial cells and the subsequent immune responses. To see the same magnitude of variation in the bacterial numbers in the model as in experiments, where the bacterial numbers range between 10 and approximately 107CFU, simulated bacterial numbers were considered in the range 10−6–1, with θBT=10−6 (refer to text 3 in the electronic supplementary material for details). For an easier comparison with the experimental growth curves, we multiplied by 107 and the logarithm of the modified is plotted in figure 2b.
The strength of phagocytosis () reveals several factors required for a correct description of the decrease in bacterial numbers. We used the Hill-type nonlinear differential equation
The maximum phagocytosis strength depends on the concentration of antigen–antibody complexes and on at the site of the infection because Ag–Ab complexes lead to better recognition of bacteria by phagocytes and IFNγ governs the rate of phagocytes' recruitment. Although IFNγ is crucial in recruiting sufficient phagocytes to the site of infection, other cytokines such as PIC contribute to a basal level, indicated here by C. V is a normalization parameter. The shape of the Hill function depends on the number of phagocytes (MP and recruited PMNs (RP)). The Hill constant H is the concentration at which the rate of increase in phagocytosis is half of the maximal rate. Thus, phagocytosis is weak when the number of phagocytes is less than the parameter H, and the sharpness of the transition between low and high phagocytosis is determined by the concentration of free bacteria multiplied by a constant k. We tuned the parameters to get the best qualitative fit of the experimental and simulated growth curves (see text 3 in the electronic supplementary material for details).
The activation of antigen-specific responses is influenced by the transport of cells and cytokines between the two compartments (see figure 1 in the electronic supplementary material). We modelled the dynamics of the transport of dendritic, Tr, Th1 and Th2 cells by Hill functions such as
(table 1). Here, the Hill constant HDC is the concentration at which the rate of transport is half of the maximal rate and the Hill coefficient n determines the slope of the curve when the rate of transport is half of the maximum. The rate of transport of cells is dependent on the interactions between surface molecules activated on the vascular walls and the cells. Hence, an initial lag due to the requirement for a high enough expression of cell-activated surface molecules can be incorporated by a sigmoidal curve with n>1.
Unlike cells, cytokines are transported to various places through the blood and lymph. Owing to the large number of unknown parameters affecting the local concentration of cytokines, we decided to randomly sample the distribution of cytokines IL4, IL12, IL10 and IFNγ over the two compartments. We introduced the random variables rc and r so that a fraction f of the total concentration of the cytokine, chosen from the interval rc−r<f<rc+r, is present in the compartment in which the cytokine is produced and a fraction 1−f is present in the other compartment. The cytokine concentrations have been experimentally observed to fluctuate. We hypothesized that these fluctuations are due to the production of cytokines by different subsets of cells in different compartments and are enhanced by the diffusion of the cytokines. To approximate the stochastic aspects of cytokine production and diffusion, we generated a new fraction f from the same interval several times during the simulation (refer to text 3 in the electronic supplementary material for the implementation of the sampling).
We modelled the differential regulation of the nodes in the two different compartments by assigning compartment-specific θ and γ values (refer to text 1 in the electronic supplementary material for a detailed description of inter-compartmental dynamics). The model generates the evolution of the states of the nodes with time, for example the time course of bacterial numbers, cytokine concentration or immune cell numbers. In §5, we compare the qualitative features of the experimental and simulated time courses.
WT and mutant infections were simulated using an initial bacterial concentration ; the initial concentrations and activities (discrete variables) of the rest of the nodes were zero. In the mutant simulation, inactive TTSS was modelled by maintaining the discrete variable of TTSS zero throughout the simulation. The concentration of each node can vary between zero and the attainable maximum values (refer to text 3 in the electronic supplementary material for a description of the attainable maximum concentrations).
The values of the parameters θ, γ, α, r, rc, n and H were selected uniformly randomly from the ranges given in table 2. We obtained these ranges by iteratively restricting the broadest relevant range (e.g. 0≤θ≤1) to improve the number of solutions. The constraints described before were used to select sets of successful parameter values. Although such sets were rare and found approximately once in 50000 simulations, we sampled the parameter space until we found 30 successful sets. The selected parameter sets outline the biologically acceptable parameter space, which we analyse in detail in §6.
Here, we describe the time course of the simulation using one representative parameter set. Certain details of the dynamics vary across different parameter sets, but all satisfy the qualitative constraints. Since the time course is the result of complex integrated behaviour of the system, the dynamics of immune components other than those specified in our constraints are predictions of the model. Note that none of the inputs to the model specify a time unit. Fixing the time unit of the model to indicate clearance at the days given in the experimental bacterial curves yielded a correspondence of the model unit time to a range of 1.3–3.5 days over the successful parameter sets. The time courses of the immune components are given until the bacteria are cleared . After bacterial clearance, most of the immune components relax into an inactive state (zero concentration) except antibodies and B cells, which remain active in all simulations, signifying immunological memory.
The model reproduced the three phases in the experimental growth curve of WT and mutant bacteria in the corresponding infections remarkably well (figure 2b). In the simulation, bacterial numbers increase exponentially, followed by a decrease in numbers due to phagocytosis. The rate of phagocytosis is dependent on the Ag–Ab complexes and IFNγ concentration and is less in phase II compared with phase III in WT infections, leading to a slower decrease in bacterial numbers in phase II. In simulation of infection by mutant bacteria, the earlier activation of IFNγ leads to a narrow phase II and faster clearance. The in vivo growth curve is a combined effect of bacterial growth and the elimination of bacteria by host immune responses. Thus, the reproduction of the growth curves indicates that the dynamic model accurately represents the host's immune response.
The experimental concentration curves of IFNγ and IL10 (see figure 2 in the electronic supplementary material) were used to set the constraints for selecting the successful parameters. Figure 3 shows the time courses of IFNγ and IL10 in WT and mutant simulations. All of the parameter sets reproduced the qualitative features formulated as constraints. In addition, in 83 per cent of the WT simulations, the concentrations of both IL10 and IFNγ are high (higher than their respective θ's) when WT bacteria are cleared. This feature of the experimental concentration curves was not explicitly included in the constraints because the regulation of IFNγ by IL10 was not clear from the experimental data. Our assumption of a cooperative inhibition of IFNγ by IL10 and the processes activated by TTSS (see text 2 in the electronic supplementary material) enables the simultaneous activity of IL10 and IFNγ in phase III. Furthermore, as observed in the experimental curves (see figure 2 in the electronic supplementary material), IL10 concentration drops after an initial increase, which is followed by another increase in the concentration level.
In the mutant simulation, the earlier production of IFNγ is by MP (node MP) and not by their other activator, Th1 cells. Indeed, several studies have speculated that the earlier production of IFNγ is not by T cells but by cells implicated in the Th1 response (Carvalho-Pinto et al. 2002; Schoenborn & Wilson 2007). We found two behaviours of the decrease in the IFNγ concentration: a drop below θIFNγI or a significant decrease that nevertheless did not fall below θIFNγI. We accepted both types of behaviour since the strength of the inhibition was not explicitly modelled. In 97 per cent of the cases, inter-compartmental dynamics were not sufficient to reproduce the drop in the IFNγ concentration. This result supports our hypothesis that the reduction in the IFNγ concentration following the peak is the outcome of a combined effect of inter-compartmental dynamics and production of inhibitory cytokines, represented by the nodes IL10 and IL4, which cooperate to inhibit IFNγ (see text 2 in the electronic supplementary material). The ratio of the day of bacterial clearance to the day on which the decrease in the concentration of IFNγ was observed is somewhere between 1.3 (28/21) and 4 (28/7) in the experimental observations, in agreement with the range of 1.1–3.3 in our simulations. Note that this ratio was not used for the parameter selection.
We used the successful parameter sets to perform perturbations/mutations experimentally studied in Pilione & Harvill (2006). Deletions were simulated by assigning zero to the discrete variable of the node throughout the simulation. The model could reproduce the higher and uninterrupted activation of IFNγ when IL10 was deleted (experimental mutant IL10−/−). We also reproduced the early clearance of bacteria in IL10 deletion (experimental mutant IL10−/−) and later clearance in IFNγ deletion (experimental mutant IFNγ−/−). Note that these time courses were not used as inputs to the model.
The time courses of immune components not directly constrained in our model are described in detail in text 4 in the electronic supplementary material. We find that several nodes can have multiple activation peaks, depending on the interplay between the half-life of the immune components and the duration of the infection. Here, we highlight some of the main features of the adaptive immune response. In agreement with the literature that Th1 cells are crucial in the clearance of micro-parasites, we see a higher Th1 response compared with Th2 response for all successful parameter sets. We find that Th2 cells are activated for a shorter time, which nevertheless is sufficient to induce B cell proliferation. IL4 is active only for a small amount of time for all parameter sets. We find that stochasticity in the distribution of cytokines over the two compartments plays an important role in the switch between the dominance of Th1- versus Th2-type immune responses. Text 5 in the electronic supplementary material shows that the switch is not possible in the absence of inter-compartmental dynamics. Saturating behaviour with different initial points was observed for the components of humoral immunity. Since the time course of most of these immune responses is not known, the simulations make relevant predictions that are corroborated by the known information.
The fact that multiple solutions with different parameter combinations exist suggests that a multitude of possible mammalian immune responses, with different rates of cytokine production or cell effector functions, can clear a Gram-negative bacterial infection of the lungs. Thus, the multiple solutions we find can illustrate the variability in the immune response across different host–pathogen systems. To investigate the existence of constraints or shared properties of these multiple solutions, we studied the existence of any restrictions in the successful parameter values (e.g. whether certain parameters need to be low) and of relationships among parameters (e.g. whether a low value of one parameter needs to be paired with a high value of another).
The θ and γ values of most of the nodes are relatively equally distributed over the range given in table 2, with a median between 0.2–0.3 and 0.4–0.6, respectively (see figure 3 in the electronic supplementary material). Figure 4a exemplifies the characteristic distribution by the nodes PIC and TrI. However, there are notable exceptions to the uniform θ distributions, e.g. the medians of nodes DCII and Th2II are lower, whereas IL4II, EC, Th1I have a higher median. The IL10I node stands out for its low median and range of lower values (figure 4a). Notable exceptions to the uniform γ distributions are DCI and PIC with lower medians, while RP and IL10I have higher medians and lower values are lacking in those nodes (figure 4b).
Both θ and γ box plots suggest the existence of subtle differential regulation for some of the parameters, either by restrictions in their range or through parameter correlations. To better understand the effect of parameters and parameter correlations on the activity of the nodes, we constructed a simple illustrative example consisting of two nodes, a source node A and a target node B. We separately considered the cases when A activates or inhibits B. We described the dynamics of the nodes by the hybrid formalism. The continuous variable of node A is described by , and node B is given by (when A activates B) or (when A inhibits B). The initial condition of A is , while in the case of activation and (i.e. the maximal value) in the case of inhibition. This example is described in detail in text 6 in the electronic supplementary material. The activity of B is dependent on the θA, θB, γA and γB parameters, as illustrated in figure 5. We first used this example to separately analyse the effect of each parameter on the activity of B to interpret the observed ranges of some θ and γ parameters. For instance, a comparison of figure 5a with 5b indicates that a higher value of θA leads to the later activation of B.
IL10 is an anti-inflammatory cytokine, which inhibits IFNγ and PIC in the model. Solution of the illustrative example (see text 6 in the electronic supplementary material) indicates that, in an inhibitory interaction, the activity of the target node is minimized for the lower θA values and higher γ values. Thus, the low range of θIL10I and high range of γIL10I indicate that IL10 effectively suppresses the activity of its targets. A similar argument applies to Th2II. Next, DCII can act as both a source and a target in activating interactions. Solution of the illustrative example indicates that the low θDCII values observed can maximize the activity of DCII's target node or the activity of DCII itself. θIL4II have higher values; IL4II can act as a source or target in inhibitory interactions and, according to the illustration, the higher threshold values imply that IL4II can be effectively suppressed, and is less effective in inhibiting its targets (corroborated by the observation that smaller concentrations of IL4 are observed in Bordetella infections). The illustration indicates that higher γ values of source and target nodes can ensure the effective activation or inhibition of target nodes. Noticeably higher ranges of γRP, which can act as both a source or target node in activating interactions, suggest the importance of the activity of RP and its targets.
To find any interdependence among parameters, we performed a correlation analysis on the successful parameter sets. A significant degree of interdependence may be a reflection of a subtle regulatory relationship not depicted in the static network or a constraint meant to ensure the reproducible activation (or suppression) of key immune processes despite variations in parameters.
We first analysed the correlations of the parameters θ and γ between different nodes and found 42 (10%) significant θ–θ correlations and 35 (8%) significant γ–γ correlations. Next, we analysed the correlations between θ of one node and γ of another to find 79 (9%) correlations. The distribution of all statistically significant correlations (p<0.05; listed in table 1 in the electronic supplementary material) among five node categories is shown in figure 4 in the electronic supplementary material. The nodes with the highest number of correlations, namely IL12I whose θ value is correlated with the γ values of 10 nodes, Cab whose θ value is correlated with the θ values of five nodes and RP whose γ value is correlated with the γ values of seven nodes, are predicted to be key regulators of the immune response.
To better understand the meaning of these correlations, we again turn to the illustrative example and study the effect of correlations on the activity of the target node B (see text 6 in the electronic supplementary material for more details). We analytically solved representative cases when γA–γB, θA–θB, γA–θB, γB–θA, γA–θA and γB–θB are positively or negatively correlated. We hypothesize that the parameter correlations are meant to optimize the effectiveness of the regulatory activity of the source node A. According to our hypothesis, a change in one parameter value may be compensated by a complementary change in the other parameter value to ensure a low value of the time point tB when the discrete variable of node B switches from its initial value. The effect of parameter correlation on the activity of B is illustrated in figure 5. In an activating interaction, an increase in θA induces an increase in tB but tB can be reduced again if θB decreases. Thus, pairing an increase in θA with a decrease in θB (in other words, a negative θA–θB correlation) may have a stabilizing effect on the activity of the target node.
In contrast to the single-parameter variations that lead to a monotonic increase or decrease in the activity of B, we find that correlated variation of two parameters, while keeping the other two constant, can lead to an extreme (maximum or minimum) activity at intermediate parameter values, and to a consistently high activity (in the case of activation) or low activity (in the case of inhibition) over a significant parameter range. Four types of correlations, namely negative θA–θB and γA–γB correlations and positive θA–γB and γA–θB correlations, maximize the activity of the target node in an activating interaction. Similarly, four types of correlations, namely positive θA–θB and θA–γB correlations and negative γA–γB and γA–θB correlations, minimize the activity in an inhibitory interaction. In the other situations, the activity is either monotonic or parameter dependent.
This simple illustration should not be extrapolated to describe the correlations between the parameters of two nodes with a complex regulatory relation. For example, several pairs of nodes with correlated parameters had both activating and inhibiting paths connecting them (e.g. θT0–θOab), and therefore the effective type of interaction between them was contingent on the time-dependent state of the nodes in these paths, which frequently included the DC, Th1, Th2, IL12 and IL4 nodes. We refer to such node pairs as mixed regulation. In many cases, paths exist in both directions (e.g. A→B and B→A); for these, we decided the effective source and target based on the directionality corresponding to the pure (not mixed) regulation or, if both regulations were pure, the regulation with the shorter path. In five cases, directionality could not be decided because the paths in both directions had equal length or a difference of 1; we denote these cases mutual regulation. We found five correlations between the parameters of TTSS and the parameters of either IL12 or IL4, although there was no effective path between TTSS and IL4 or IL12. These correlations suggest that the existence of subtle interdependencies was not incorporated in the interaction network. For example, both TTSS and IL4 activate IL10, which inhibits PIC, whereas through another pathway TTSS activates PIC; thus the correlation among the TTSS and IL4 parameters suggests a complex regulation of host processes by TTSS.
Overall, 43 per cent of the statistically significant source–target correlations agree with the illustrative example (refer to text 6 in the electronic supplementary material for the description of θ–θ, γ–γ and γ–θ correlations), and none of the remaining cases contradict it, lending support to the hypothesis that the role of parameter correlations is to ensure the effectiveness of the regulatory interaction between the source and target node. For several significant correlations that do not agree with our illustration, the two-node premise of the illustration was not a good fit owing to combinatorial regulation of the target node. For example, the rule for BC is BC=Th2 or BC, thus Th2 is required only for first activation. Thirty-nine parameter correlations corresponded to node pairs with mixed regulation, including 16 correlations among the nodes DC, Th1, Th2, IL12 and IL4 that are connected by both mutual and mixed relationships. The scarcity of γ–γ correlations among these nodes suggests that the activation threshold θ plays an important role in the decision of the differentiation of T helper cells. Some of the correlations among mixed regulated nodes indicate an interesting interdependence; for example, γEC–θIL12I suggests the regulation of Th1 responses by innate immune responses. Compartment-specific node correlations were also observed (e.g. θDCI is correlated with θIL12II but not with θIL12I).
α is a scaling factor that governs the rate of activation and decay of cytokines. We found that the α values of all cytokines have a median of 2. αPIC and αIFNγ have higher values than the other α's (see figure 3c in the electronic supplementary material), indicating that PIC and IFNγ tend to be more variable than the other cytokines. This is consistent with the inhibitory functions of IL10, IL4 and IL12. We found a significant degree of positive correlation among the α parameters of cytokines in the T helper cell network, while PIC is correlated (negatively) only with IL12.
The transport of cells between the two compartments is described by the H and n parameters, while the r and rc parameters describe the inter-compartmental dynamics of four cytokines. A larger range of HTh1 than that of HTh2 (figure 6a) indicates a higher variation in the concentration required for effective transport of Th1 cells. Correlations between the H, n, r and rc of various nodes suggest the existence of interdependences among the flow of cells or cytokines between the two compartments. The two significant H correlations were of HTr with HTh1 and HTh2, both negative correlations, suggesting that T regulatory cells modulate the transport of both Th1 and Th2 cells. The box plot in figure 6b indicates that Th2 and DC frequently have smaller Hill coefficients n, whereas Tr tends to have higher values, indicating that Tr transport is more switch-like than the others. The only significant correlation was a negative correlation between nTh1 and nTh2, supporting the conclusion that inter-compartmental dynamics, in addition to mutual inhibitory interactions, play a role in the Th1–Th2 dichotomy. The positive (HTr–nTh2) and negative (HTh1–nTr) correlations between the H and n parameters indicate that the nonlinearity (determined by n) of the rate of flow from one compartment to the other of Th2 and Tr is influenced by the concentration of Tr and Th1, respectively. Moreover, the H and n parameters (HTr–nTr) of Tr also influence each other.
In the model, on average, a fraction rc is present in the compartment in which the cytokine is produced, a fraction 1−rc is present in the other compartment and the variability in the concentration is proportional to r. The box plot in figure 6c shows that the rc values vary more among nodes than the r values. The median of rcIL4 is smaller than 0.5, whereas the medians of rcIFNγ, rcIL12, rcIL10 are larger than 0.5. In addition to negative correlations between rc and r of the same cytokine (expected since r<rc), we found only one positive correlation, rcIL4–rIL4. As the range of rcIL4 is below 0.5, the positive rc–r correlation of IL4 indicates a more stable (less variable) activation of IL4. Indeed, IL4 is known to be activated at low but steady amounts. We also observed negative correlation between the same parameter of different cytokines, rIFNγ–rcIL12 and rcIL10–rcIL12, suggesting that they affect each other's concentrations in the compartment.
Thus, in many instances, the correlation analysis provided results that agree with our network topology. In other cases, the analysis suggests novel regulatory relationships that provide directions for future experiments. Some of the predicted regulatory relationships agree with previously proposed relationships that nevertheless did not have enough support in the Bordetella-related literature to be included in our network, or appear as one path in a mixed regulation. Since the mixed regulations describe context-dependent interactions, these correlations suggest that more detailed experiments are needed to reveal the full extent of the regulation. For example, we found many correlations between the parameters of epithelial cells and Th2-related components (e.g. θEC–θTh2I, γTh2I–θEC, γEC–θIL4II). It has been shown by independent studies that cytokines produced by epithelial cells in response to bacterial pathogens, e.g. IL6 (Agace et al. 1993; Weinstein et al. 1997), can promote Th2 differentiation and inhibit Th1 differentiation (Diehl & Rincon 2002). We also found correlations between IL12 and Cab (θCab–θIL12I, γCab–θIL12I), the implication of which, that Th1-related cytokines promote complement-fixing antibodies, has been proposed (Janeway et al. 2001; Peng et al. 2002). Furthermore, the role of dendritic cells and T0 cells in differentially regulating the Th1 or Th2 response has been observed in many infections (Kalinski et al. 1999; Reis e Sousa et al. 1999; MacDonald & Pearce 2002), and is also reflected in the correlations.
The immune response can be thought of as a dynamic interplay among interdependent interaction networks such as transcription factor–gene interaction networks, cell–cell interactions and cytokine interaction networks. Traditionally, models of the immune system focus on subsystems composed of less than five components. Although properly chosen subsystem models have been successful, they cannot address the emergence of holistic properties. Instead, it may be advantageous to coarse-grain the components (network nodes) such that what is represented as the response of one node is in reality a collective response of some of the networks mentioned previously. The resulting network model, such as the one presented here, not only provides holistic insight but can also predict the qualitative outcome of the subsystems. Since the input to the network assembly process is from completely independent studies, the fact that the dynamic model of the network reproduces experimental observations provides verification for the network and the dynamic model. The evaluation of the data in the form of qualitative constraints allowed us to overlook the variations across experiments and laboratories. Moreover, this model enables us to extend the results based on an experimental system to different host species and pathogen strains. Additionally, the model separately incorporates local and systemic immune responses to pathogens, an issue which is rarely discussed in immunological studies.
Experimental results often infer the regulation between two components, e.g. IL10 and IFNγ, at a certain time point. However, this regulation may be context dependent, and therefore iterative network improvements towards the topology that describes the observed dynamic outcome may reveal the existence of co-regulators. Indeed, our extended network indicates that IL10 is not the only cytokine responsible for inhibiting IFNγ but other Th2-related cytokines such as IL4 may also play a regulatory role in this inhibition. This is supported by the observation that IL10 and IFNγ can be active at the same time. The optimal topology of the sub-network of T helper cell differentiation suggests that IL12 is necessary for the differentiation of T0 cells into Th1 cells and its absence leads to differentiation into Th2 cells. This is supported by the fact that IL4 is produced by differentiated Th2 cells. Although this topology describes the known observations, cell surface molecules on antigen-presenting cells and their interactions with specific antigens also play a role in directing T helper cell responses (Zhang et al. 2004). The network also suggests several interactions of TTSS with other components such as the induction of IL10 by splenocytes (MP). These model predictions can serve as targets for future experiments.
Most theoretical models focus on experimentally well-studied cellular pathways that have relatively few components, and they use ODEs describing the kinetics for the production and decay of all components. The drawback of the differential equation-based method is the necessity to know or estimate the kinetic parameters of the underlying processes. These parameters are often not available from immunological experiments. On the other hand, qualitative data, which are usually more abundant, are underestimated in their ability to provide mechanistic knowledge and constructive hypotheses. The evidence that the input–output curves of many regulatory relationships are strongly sigmoidal and the evidence that regulatory networks maintain their function even when faced with fluctuations in components and reaction rates (von Dassow et al. 2000; Alon 2006) support the applicability of qualitative models. The success of these models in describing specific cellular systems and processes such as flower development (Espinosa-Soto et al. 2004), the yeast cell cycle (Li et al. 2004) and Drosophila embryonic development (Sanchez & Thieffry 2001) indicates that the topology of regulatory networks has a significant role in restricting their dynamical behaviour.
In our previous study, we developed a Boolean model, which reproduced the clearance or persistence of bacteria under various conditions and indicated the importance of the degradation of cytokines and of bacterial virulence factors, and the existence of activation thresholds. In the present study, the activation threshold θ and decay rate γ, respectively, make predictions about the relative concentrations required for the activation of downstream processes and about the decay time scales. We use node-specific values for these parameters, which are more detailed than the traditional assumption of θ=0.5 and γ=1 or the recently introduced variation of scaling the time units in each differential equation (Chaves et al. 2006).
To determine whether the model parameters need to be fine-tuned, we used a random sampling method. Although the solutions were rare, we found several of them, and the values of most of the parameters were evenly distributed over their respective ranges. This result corroborates previous conclusions that the topology of regulatory networks is more important than the fine-tuning of kinetic parameters (von Dassow et al. 2000). It will be interesting to explore whether the incorporation of more experimental data, as they become available, will narrow down the successful parameter space. In accordance with previous observations (Chaves et al. 2006), we found that the optimal threshold for turning a node's activity on is between 0 and 50 per cent of its maximal concentration. Unlike other nodes, the lower bound of the θIL12II and θIL4II range was 6 per cent and 11 per cent of the maximal concentration, respectively. These results indicate that higher amounts of IL12 and IL4 need to be produced in the lymphatic tissues where Th cell differentiation takes place.
The predicted time courses of our model confirm and provide further detail for the time courses of our previous discrete model (refer to text 7 in the electronic supplementary material for a detailed comparison). The discrete model represented the outcome of the mutual inhibition between Th1 and Th2 cells (and related processes) as a bistable switch and revealed that decay is required for the switching behaviour. Deeper analysis by the hybrid model shows that the dominance of one T helper cell type over the other depends on the local cytokine environment. The model indicates the importance of inter-compartmental dynamics. We observed different ranges of parameters in the two compartments and interestingly found strong correlation between the rc parameters of IL12 and IL10. The model also suggests the existence of differential regulation of Th1 and Th2 cells by other immune components through the existence of differential correlations, e.g. θDCI is correlated with but not with . Thus, the current model indicates that several factors, including various immune components and inter-compartmental dynamics, play a role in modulating T helper cell response.
The logical network representation and intuitive dynamical model will allow microbiologists and immunologists to see the knowledge gaps in the results from different laboratories and appreciate the synergy between considering both the pathogen and the host. The compartmentalization of the model allows its extension to study co-infections in which two pathogens infect different organs of the same host. Since the hosts studied in laboratory experiments are not always the natural host of the pathogen, models such as this are crucial to explore the outcomes of the disease in other hosts. Bordetella bronchiseptica is the ancestor of two human respiratory pathogens, B. pertussis and B. parapertussis. Our future goals include the development of a hybrid model for B. pertussis. Key requirements for this model are the availability of comparative input data (similar to the IFNγ time courses in WT and TTSS deficient B. bronchiseptica), and of evidence of the mechanisms by which B. pertussis is modulating the host's immune system. Owing to the complex regulation of immune components by the virulence factors of B. pertussis, establishing causal relationships between virulence factors and immune components is a challenge and future targeted experiments will be required to do so. In general, species-specific factors can be readily incorporated in the network, allowing the comparison between different species to recognize crucial virulence mechanisms and identifying targets for intervention.
This work was partially supported by NIH grant 1-R56-AI065507-01A2 to E.T.H., NSF grant CCF-0643529 to R.A. and an Irvington Institute and Cancer Research Institute fellowship to J.T. We are grateful to István Albert for developing the software BooleanNet that forms the basis of our code.
Note that this meaning of ‘phases’ is different from its usage in the expressions ‘Bvg− phase’ (non-infectious) or ‘Bvg+ phase’ (infectious) of bordetellae.
Supplementary texts and figures are included in one file