Gastrointestinal helminths increase Bordetella bronchiseptica shedding and host variation in supershedding

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

218 views

Version of Record published: November 8, 2022 (This version)
Accepted: October 11, 2022
Received: May 14, 2021
Preprint posted: May 6, 2021 (Go to version)

1. Of interest
Multiple preferred escape trajectories are explained by a geometric model incorporating prey’s turn and predator attack endpoint

Yuuki Kawabata, Hideyuki Akada ... Paolo Domenici

Research Article Updated Mar 31, 2023
Further reading

Article
Figures and data
Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

Co-infected hosts, individuals that carry more than one infectious agent at any one time, have been suggested to facilitate pathogen transmission, including the emergence of supershedding events. However, how the host immune response mediates the interactions between co-infecting pathogens and how these affect the dynamics of shedding remains largely unclear. We used laboratory experiments and a modeling approach to examine temporal changes in the shedding of the respiratory bacterium Bordetella bronchiseptica in rabbits with one or two gastrointestinal helminth species. Experimental data showed that rabbits co-infected with one or both helminths shed significantly more B. bronchiseptica, by direct contact with an agar petri dish, than rabbits with bacteria alone. Co-infected hosts generated supershedding events of higher intensity and more frequently than hosts with no helminths. To explain this variation in shedding an infection-immune model was developed and fitted to rabbits of each group. Simulations suggested that differences in the magnitude and duration of shedding could be explained by the effect of the two helminths on the relative contribution of neutrophils and specific IgA and IgG to B. bronchiseptica neutralization in the respiratory tract. However, the interactions between infection and immune response at the scale of analysis that we used could not capture the rapid variation in the intensity of shedding of every rabbit. We suggest that fast and local changes at the level of respiratory tissue probably played a more important role. This study indicates that co-infected hosts are important source of variation in shedding, and provides a quantitative explanation into the role of helminths to the dynamics of respiratory bacterial infections.

Editor's evaluation

The authors perform experimental infections with rabbits to study how coinfection with one or more helminths affects the shedding of the respiratory bacterium Bordetella bronchiseptica. The results show that shedding varies strongly from one individual to the next and that co-infections with helminths lead to increased levels of shedding. The authors nicely combine within-host kinetics modelling and their longitudinal data to estimate key parameter values associated with bacterium and immune growth rates in the four conditions. These suggest that the shedding differences can be explained by differences in bacterial growth.

https://doi.org/10.7554/eLife.70347.sa0

Introduction

Individual variation in pathogen transmission can increase the basic reproduction number R₀ of a pathogen and determine whether an infection will invade and spread or stutter and quickly fade out in a population of susceptible hosts (Keeling et al., 2001; Lloyd-Smith et al., 2005). One of the causes of this variation is associated with differences in the amount and duration of pathogen shedding, whereby some infected hosts shed disproportionately more and for longer than the average population, the so called supershedders (Chase-Topping et al., 2008; Gopinath et al., 2013), while others do not shed at all (Chen et al., 2006; Hadinoto et al., 2009; Leung et al., 2015). Co-infected hosts are frequently proposed to contribute to this variation, which could emerge from interactions between pathogen species mediated by the host immune response and the consequences on host infectiousness (Sheth et al., 2006; Graham et al., 2007; Richard et al., 2014).

While studies on the immunological response to multi-species infections has provided insight to the interactions between the host and its pathogens, there remains a need to identify how these processes relate to onward transmission, specifically the patterns of bacterial shedding. Dynamical mathematical models are particularly useful in disentangling these complexities as they can generate mechanistic-driven hypotheses that can be examined in relation to empirical data (Smith et al., 2013; Byrne et al., 2019). In the current study, we applied this general approach to investigate the dynamics of shedding of the respiratory bacterium Bordetella bronchiseptica in rabbits experimentally co-infected with one or both of the gastrointestinal helminths Trichostrongylus retortaeformis and Graphidium strigosum. We explored to what extent helminths could alter the level of B. bronchiseptica shedding over time, whether the trend varied depending on the helminth species and to what extent the host immune response could explain the patterns observed.

B. bronchiseptica is a highly contagious bacterium of the respiratory tract that causes multiple symptoms and infects a wide range of mammals (Goodnow, 1980). In rabbits and mice, and likely other mammal species, B. bronchiseptica is removed from the lower respiratory tract (lungs and trachea) via phagocytosis stimulated by a Th1 inflammatory reaction that involves cell mediated antibodies and neutrophils (Thakar et al., 2007; Thakar et al., 2009; Thakar et al., 2012). Bacteria persist in the nasal cavity, although they are partially reduced by IgA antibodies in naïve and immunized mice (Kirimanjeswara et al., 2003; Wolfe et al., 2007). Transmission is poor among wild-type laboratory mice but increases among TLR4-deficient mice (Rolin et al., 2014). The TLR4-deficient response is associated with neutrophil infiltration, and the intensity of shedding has been found to be positively correlated with neutrophil counts. In contrast, rabbits naturally shed (Pathak et al., 2010) and efficiently transmit B. bronchiseptica between animals (Long et al., 2010). The evidence of rapid outbreaks in pet kennels, livestock holdings and laboratory rabbitries are consistent with rapid transmission among animals in close contact. Occasionally, B. bronchiseptica spills-over into humans but there is no evidence of sustained onward transmission, as these human infections have invariably occurred in immunocompromised individuals (Goodnow, 1980; Woolfrey and Moody, 1991). Humans are primarily infected by the species B. pertussis and B. parapertussis which are responsible for regular whooping cough outbreaks worldwide (Domenech de Cellès et al., 2018). Given the close relatedness between Bordetella species, and considering the many similarities in the kinetics of infection and the immune response, B. bronchiseptica provides a good model system to explore the dynamics of shedding in Bordetella infections and the interaction with other infectious agents, like gastrointestinal helminths.

In co-infections with other respiratory pathogens B. bronchiseptica contributes to exacerbate respiratory symptoms, including the development of acute pulmonary disease and bronchopneumonia, and ultimately host death (Brockmeier et al., 2008; Loving et al., 2010; Schulz et al., 2014; Kureljušić et al., 2016; Hughes et al., 2018). Parasitic helminths commonly stimulate a Th2 anti-inflammatory immune response that interferes with the Th1 response developed against B. bronchiseptica and related species (Brady et al., 1999; Pathak et al., 2012; Thakar et al., 2012). We recently showed that rabbits co-infected with either T. retortaeformis or G. strigosum carried higher B. bronchiseptica infections in the nasal cavity but there were no significant differences in the size and duration of infection in the lungs, when compared to rabbits infected with bacteria alone (Pathak et al., 2010; Pathak et al., 2012; Thakar et al., 2012). These two helminths stimulate a similar immune response, however, while the former is reduced or cleared from the small intestine, the latter persists with high intensities in the stomach (Cattadori et al., 2005; Cattadori et al., 2008; Cattadori et al., 2019). The modeling of the immune network in B. bronchiseptica-T. retortaeformis rabbits and B. bronchiseptica alone animals suggested that neutrophils and antibodies IgA and IgG are important to bacterial clearance from the lungs (Thakar et al., 2012). Most likely, this mechanism could also explain patterns observed in B. bronchiseptica-G. strigosum infections where bacteria are removed from the lower but not the upper respiratory tract (Pathak et al., 2012). Laboratory infections of mice with B. bronchiseptica supported the role of these three immune variables to bacterial clearance from lungs and trachea, and reduction in the nasal cavity (Kirimanjeswara et al., 2003; Wolfe et al., 2007).

To investigate the impact of these two helminths to B. bronchiseptica shedding, four types of laboratory infections were investigated: i- B. bronchiseptica (B) only, ii- B. bronchiseptica-G. strigosum (BG), iii- B. bronchiseptica-T. retortaeformis (BT), and iv- B. bronchiseptica-T. retortaeformis-G. strigosum (BTG). For the first three experiments, studies on the dynamics of infection and related immune responses are available in Pathak et al., 2010; Pathak et al., 2012; Thakar et al., 2012 and key findings have been described in the previous section. The triple infection (BTG) followed the same experimental design and laboratory procedures, and is presented here for the immune components relevant to this study, together with novel data on B. bronchiseptica shedding from all the four types of infection (details in Materials and methods). Longitudinal data on immunity and bacteria shedding over four months were then used to develop an individual-based Bayesian model that was applied to each type of infection, independently. First, we built a dynamical model that investigated the relative contribution of neutrophils and specific IgA and IgG to bacterial infection in the whole respiratory tract. Second, we applied this dynamical model to laboratory data, where shedding was linked to the estimated bacterial infection. Helminths were assumed to alter the immune response to B. bronchiseptica by impacting the magnitude and time course of the three immune variables. Simulations provided a parsimonious explanation of the dynamics of shedding and how the two helminth species, with contrasting dynamics, contributed to variation between types of infection and between individuals within the same infection.

Results

Experimental results: Helminth infections boost B. bronchiseptica shedding

For each type of infection, the number of bacteria shed, as determined by contact with a BG-blood agar petri dish for a fixed time (CFU/s), were collected every week or two-three times a week for every rabbit (see Materials and methods). There was large variation in the level at which B. bronchiseptica was shed both within and between the four types of infection (Figure 1). Some rabbits consistently shed large amounts of CFUs while others were low shedders or did not shed at all, specifically: 3 out of 16 for B. bronchiseptica only animals (B), 1 out of 23 for B. bronchiseptica-G. strigosum (BG), 3 out of 20 for B. bronchiseptica-T. retortaeformis (BT) and 0 out of 24 for the triple infection (BTG). Since we were interested in the dynamics of shedding, these non-shedders were excluded from the subsequent analyses. High shedders were found during the first 60 days post infection, although some of these rabbits continued to shed a large number of bacteria well beyond this time (Figures 1 and 2). The median trend of each group reflects the large number of events at medium-low intensity and the large variation among individuals observed. The level of shedding was significantly higher in rabbits co-infected with helminths than B. bronchiseptica only animals (mean shedding (CFU/s) and 95% CI, BTG: 0.42, 0.31–0.54; BG: 0.32, 0.25–0.41; BT: 0.26, 0.18–0.35 and B: 0.02,–0.01-0.03; Kruskal-Wallis test with Bonferroni correction for multiple testing, for all: p<0.05). When we considered the co-infected hosts, no significant differences were observed among the three groups (p>0.05) although the triple infection showed a tendency for higher shedding (Figure 2). Across the four infections up to 55% of the shedding events were classified as null, where rabbits interacted with the petri dish but did not shed even though they were known to be infected (Figure 3).The highest percentage was found in the B. bronchiseptica only group (56.2%), followed by the dual infections (BT: 40.3% and BG: 29.8%) and last the triple infection BTG (25.8%) (Figure 3). This latter group showed the largest variation in the level of shedding (Negative Binomial k and 95% CI, BTG: 0.21, 0.16–0.26) followed by the dual infections BG (0.37, 0.29–0.45) and BT (0.45, 0.33–0.57); for all, the distribution was not significantly different from the frequency data (p>0.05). There were only three bins to represent the frequency of shedding in the B. bronchiseptica alone group and the Negative Binomial distribution could not be fitted. These findings suggest that the presence of helminths increased the level of bacteria shed in the environment and the frequency of these events. Having helminths also increased the variation in the magnitude of these events among rabbits within the same group as well as between groups. Together these findings indicate that co-infections with helminths, irrespective of the species, should lead to a higher probability of B. bronchiseptica onward transmission given a contact.

Figure 1

Download asset Open asset

Observed levels of B. *bronchiseptica* shedding by contact of infected rabbits with a Bordet-Gengou agar petri dish.

Rabbits (infected shedders and non-shedders) from the four types of infection are reported, *B. bronchiseptica* (B), *B. bronchiseptica-G. strigosum* (BG), *B. bronchiseptica-T. retortaeformis* (BT) and *B. bronchiseptica-T. retortaeformis-G. strigosum* (BTG) (top-down). Each infection is presented following the experimental design where rabbits are ordered by their time of sacrifice (early/left, late/right) during the infection (week 4^th to 21^st, dotted vertical red lines). For every rabbit the following are reported: observed shedding events (black points), 25^th and 75^th percentiles (top and bottom hinges), median values (thick horizontal lines) and values within the 1.5 inter-quartile range (IQR, whiskers extending up to 1.5*IQR). To facilitate the visualization, B infection is reported on a smaller y-axis. The level of shedding is quantified as number of Colony Forming Units per second (CFU/s).

Figure 2

Download asset Open asset

Observed level of B.bronchiseptica shedding by days post infection (Dpi) and type of infection.

Shedding events (points), host individual trajectories (thin lines) and median population trend (smoothed thick lines) with 95% CI (blue shadows) are reported. The trajectories of few individuals are highlighted with different colors to facilitate visualization. The abrupt interruption of individual trends is caused by the removal of rabbits at fixed time points. To reduce overcrowding, host trajectories with at least four shedding points larger than zero are included. Note the much lower shedding values, and different y-axis, for the B group. Further details in Figure 1.

Figure 3

Download asset Open asset

Frequency of observed B.bronchiseptica shedding events for the four types of infection.

Events have been grouped by shedding level using 0.5 CFU/s unit intervals, the first interval represents events with zero shedding. Red lines indicate the 95^th and 99^th percentile thresholds, respectively, estimated using the entire experimental dataset; cases above the threshold value represent supershedding events. A Negative Binomial distribution is fitted to each co-infection group (blue line), while it was not possible for the B group due to the few bins.

Experimental results: Helminth infections promote B. bronchiseptica supershedding

Few events were characterized by high bursts, very large numbers of CFUs counted in a petri dish (Figures 2 and 4). The definition of supershedder, namely, the threshold above which there is evidence of a supershedding event, depends on the pathogen and the host, and it is usually based on the assumption that hosts carry one single infection (Chase-Topping et al., 2008; Gopinath et al., 2012). If we estimate this threshold using B. bronchiseptica only infected animals, and define supershedders as the hosts that have at least one shedding event above this threshold, for example the 99th percentile, our cut-off value is CFU/s=0.27 (Figure 3). This indicates that none of the rabbits with B. bronchiseptica (0 out of 13) can be classified as supershedders. However, if we apply the same cut-off to the co-infected groups the percentage is considerably high, namely: 59% (13 out of 22) for BG, 53% (9 out of 17) for BT and 42% (10 out of 24) for the BTG group. This threshold definition provides a common reference value for our infections, however, it is dependent on host status (i.e. calculated using hosts with single infections), and is not representative in other settings, such that it overestimates the number of supershedders in rabbits that are co-infected. A way to overcome some of the limitations of this approach is to estimate the cut-off value using the whole dataset, irrespective of their type of infection. When we apply this approach the 99th percentile threshold is now CFU/s=8.58 and the fraction of hosts with, at least, one supershedding event above this value is: 0% (0 out of 13) for B, 32% (7 out of 22) for BG, 18% (3 out of 17) BT and 21% (5 out of 24) for BTG (Figure 3). Similarly, if we select the less limiting 95th percentile threshold, the new common value decreases to CFU/s=2.71 and the number of rabbits with at least one supershedding event is still zero for B but increases to 45% (10 out of 22) for BG, 29% (5 out of 17) for BT and 29% (7 out of 24) for BTG (Figure 3). Consistent with the general pattern of shedding, supershedding events were found more often in the initial four/five weeks post infection (Figures 1 and 2), although later events were also observed, especially for the BG and BTG groups, and with some rabbits contributing multiple times. Importantly, since co-infected rabbits have higher levels of shedding and a much higher probability of becoming supershedders (whether this is at the 95th or 99th percentile threshold), we should expect a higher probability of onward transmission given a contact during these disproportionate events, compared to B. bronchiseptica only rabbits or rabbits that have average values of shedding.

Figure 4

Download asset Open asset

B. *bronchiseptica* shedding on BG-blood agar petri dishes.

Examples of A: supershedding event and B: average shedding event.

Modeling results: Changes in the relative role of neutrophils and antibodies contribute to variation in B. bronchiseptica dynamics of infection between groups

We examined the hypothesis that B. bronchiseptica shedding is controlled by neutrophils, species-specific IgA and IgG produced against the bacterial infection in the respiratory tract, and helminths affect shedding by altering the magnitude and time course of the three immune variables. Our model formulation did not explicitly include the intensity of infection of the two helminths, rather, we examined how they impacted the immune response and its consequences on shedding. Here, we report on the relationships between immune response and dynamics of infection from the dynamical model fitted to the experimental individual data, while in the next section we examine the relationship between predicted level of infection and estimated intensity of bacteria shedding. To provide meaningful results, and to avoid issues associated with model convergence, individual time series with three or less shedding events larger than zero were omitted from the modeling. Similarly, for B. bronchiseptica only rabbits model fitting was performed on the data pooled at the group level because of the small number of hosts with enough shedding events larger than zero (details in Materials and methods). The description of the model framework, fitting and validation, including parameter calibration, model selection and sensitivity analysis of key parameters, are outlined in Materials and methods and Appendix. Consistently among the four types of infection, simulations indicated that there was a rapid bacterial replication following the initial inoculum, the growth rate, r, was high for BTG (group level value and [95% CI]: 1.47 [1.10–1.83]) and BG (1.18 [0.83–1.54]), less so for BT (1.01 [0.83–1.18]) and much lower for B single infection (0.06 [0.006–0.12]) (Tables 1 and 2). The growth of B. bronchiseptica prompted a rapid immune response. Model fitting well captured the empirical trends of neutrophils and specific antibodies over time, both as a profile representative of each type of infection and as a trend of every rabbit within each group (Table 2, Figure 5). In all four groups, neutrophils showed the fastest rate of increase, a₁, with a tendency to be faster in B. bronchiseptica only rabbits, followed by specific IgA and then specific IgG (Table 2, Figure 6). Neutrophils were also the fastest to decrease, at a rate b₁, and to return to original values with the neutralization, albeit no complete clearance, of bacteria from the respiratory tract (Figures 5 and 6). As reported in our previous studies of this system, B. bronchiseptica is removed from lungs and trachea but persists in the nasal cavity (Pathak et al., 2010; Pathak et al., 2012; Thakar et al., 2012), a pattern also observed in mice (Kirimanjeswara et al., 2003). The increase of specific IgA and IgG was slow, however, once IgG reached the asymptotic value at around 20 days post infection, it remained high throughout the trials, while IgA peaked around day 20 and slowly decreased thereafter, although it had not fallen to baseline levels by the end of the trials (Figures 5 and 6). The relative contribution of the three immune variables to B. bronchiseptica neutralization, c₁, c₂, c₃, varied among the four types of infection. Model fittings suggested that neutrophils primarily controlled bacteria in the single infection, whereas a combination of neutrophils and specific IgG were probably important in BTG rabbits (Table 2, Figure 6). The coupling of IgG and IgA could contribute to explain the neutralization of bacteria in the BT group while neutrophils and IgG appeared to be relevant in the BG group. We note that timing of these interactions was also critical. Given the early reaction of neutrophils, they are expected to primarily control B. bronchiseptica early on in the infection, while specific IgG and, secondly, specific IgA appear to be more important at a later time, as they decreased more slowly, such is the case for IgA, or remained consistently high over the course of the infection, like for IgG (Figures 5 and 6). An investigation of the posterior parameter estimations confirmed that the immune response changed both between individuals within a group and between groups (Appendix 1).

Figure 5

Download asset Open asset

Estimated immune responses in the blood by time and type of infection.

Empirical longitudinal data (points) with estimated trajectories (black lines) for every individual rabbit, and group means (blue lines) with 95% CIs (blue shadows) are reported. The abrupt interruption of individual data is caused by the removal of rabbits at fixed time points. For the B rabbits we only report the trend at the group level.

Figure 6

Download asset Open asset

Estimated rates of infection in the whole respiratory tract and immune response in the blood by time and type of infection.

Means (continuous lines) and 95% CIs (shaded areas) are reported for neutrophils (black), IgA (blue), IgG (red), and *B. bronchiseptica* infection (green). Individual trends (dotted lines) have been included for the neutralization rates of the co-infected groups but not for the other rates, as CIs are very narrow and individual trends are difficult to disentangle. The abrupt interruption of individual trends are caused by the removal of rabbits at fixed time points. Full details on the estimated rates, along with their credible intervals, are available in Table 2.

Table 1

Description and unit of parameters and variables used in the dynamical and observed models.

Parameters/	Descriptions	Unit	Value
Variables
Rabbit i^th specific parameters
$a_{1, i}, b_{1, i}, c_{1, i}$	Growth, decay and B. bronchiseptica neutralization rates of neutrophils	day^-1	Estimated
$a_{2, i}, b_{2}, c_{2, i}$	Growth, decay and B. bronchiseptica neutralization rates of IgA	day^-1	Estimated
$a_{3, i}, b_{3}, c_{3, i}$	Growth, decay and B. bronchiseptica neutralization rates of IgG	day^-1	Estimated
r_i	Per capita B. bronchiseptica replication rate	day^-1	Estimated
$x_{1, i} (t)$	Time-dependent neutrophil response	cells/ml	Calculated
$x_{2, i} (t)$	Time-dependent IgA response	O.D. index	Calculated
$x_{3, i} (t)$	Time-dependent IgG response	O.D. index	Calculated
$y_{i} (t)$	Time-dependent B. bronchiseptica infection intensity	CFU/g count	Calculated
Group parameters
$μ_{a_{1}}, μ_{b_{1}}, μ_{c_{1}}$	Mean of neutrophil growth, decay and neutralization rates	day^-1	Estimated
$σ_{a_{1}}, σ_{b_{1}}, σ_{c_{1}}$	Standard deviation (S.D.) of neutrophil growth, decay and neutralization rates	day^-1	Estimated
$μ_{a_{2}}, μ_{c_{2}}$	Mean of IgA growth, and neutralization rates	day^-1	Estimated
$σ_{a_{2}}, σ_{c_{2}}$	Standard deviation (S.D.) of IgA growth, and neutralization rates	day^-1	Estimated
$μ_{a_{3}}, μ_{c_{3}}$	Mean of IgG growth, and neutralization rates	day^-1	Estimated
$σ_{a_{3}}, σ_{c_{3}}$	Standard deviation (S.D.) of IgG growth, and neutralization rates	day^-1	Estimated
$μ_{r}, σ_{r}$	Mean and S.D. of B. bronchiseptica replication rates	day^-1	Estimated
$x_{1}^{*}$	Neutrophils level at equilibrium	cells/ml	Observed
$μ_{y_{0}}$	B. bronchiseptica infection intensity at time t=0	CFU/g count	Estimated

Table 2

Model prior and posterior parameters for the immune response and B. bronchiseptica infection.

We report: mean and S.D. of priors, and mean, 95% CI and S.D. from posterior distributions for each group. The complete description of the parameters is provided in Table 1.

Parameters	Priors	Posteriors B	Posteriors BG	Posteriors BT	Posteriors BTG
$μ_{a_{1}} * 10^{- 4}$	0.10 (0.05)	2.89 (1.17–4.61)	1.08 (0.31–1.85)	3.32 (1.73–4.90)	0.93 (0.35–1.50)
$μ_{a_{2}} * 10^{- 5}$	1.10 (0.20)	2.00 (1.65–2.35)	0.81 (0.51–1.11)	0.99 (0.61–1.36)	1.36 (0.91–1.81)
$μ_{a_{3}} * 10^{- 5}$	0.30 (0.10)	1.51 (1.27–1.75)	0.66 (0.43–0.89)	0.78 (0.49–1.36)	0.86 (0.59–1.81)
$μ_{b_{1}}$	0.50 (0.10)	1.20 (0.83–1.57)	0.54 (0.37–0.71)	0.86 (0.61–1.11)	0.58 (0.40–0.77)
$b_{2} * 10^{- 3}$	0.30 (0.10)	7.07 (5.55–8.59)	3.10 (2.4–3.80)	5.57 (4.97–6.16)	7.43 (6.50–8.36)
$b_{3} * 10^{- 4}$	0.10 (0.07)	0.93 (0.03–1.83)	0.66 (0.02–1.30)	0.39 (0.01–0.77)	0.49 (0.01–0.97)
$μ_{c_{1}}$	0.20 (0.05)	0.11 (0.06–0.16)	0.13 (0.06–0.20)	0.05 (0.03–0.07)	0.30 (0.19–0.42)
$μ_{c_{2}}$	0.12 (0.02)	0.044 (0.003–0.085)	0.13 (0.03–0.22)	0.22 (0.16–0.27)	0.01 (.0004–0.02)
$μ_{c_{3}}$	0.10 (0.07)	0.042 (0.002–0.083)	0.29 (0.18–0.40)	0.33 (0.25–0.42)	0.32 (0.23–0.42)
$r$	1.70 (0.30)	0.06 (0.01–0.12)	1.18 (0.83–1.54)	1.01 (0.83–1.18)	1.47 (1.10–1.83)
$μ_{y_{0}} * 10^{3}$		0.68 (0.68–0.68)	0.71 (0.35–1.07)	0.68 (0.36–1.01)	0.82 (0.46–1.18)
$x^{*}$		0.05	1.89	1.93	2.75
$σ_{x_{1}}$		0.45	0.38	0.53	0.70
$σ_{x_{2}}$		0.22	0.40	0.33	0.22
$σ_{x_{3}}$		0.08	0.10	0.15	0.15
$σ_{y}$		1.00	6.02	4.50	4.01

Modeling results: immune mediated B. bronchiseptica infection explains variation in the dynamics of shedding between groups

For each individual rabbit, we examined the dynamics of shedding by linking the empirical shedding of every individual to its estimated intensity of infection in the respiratory tract through the zero inflated log-normal function and, from this, generated the simulated individual shedding data (details in Materials and methods). As previously noted, model fitting was performed on the whole dataset for the B. bronchiseptica only rabbits. Simulations showed that the rapid growth of B. bronchiseptica in the respiratory tract led to the rapid shedding in the environment (Figure 7). The estimated time to reach the peak of shedding was between 10 and 12 days post infection for the co-infected groups, specifically: 9 days for BTG, 11 days for BT and 12 days for BG. For B. bronchiseptica only rabbits, technical problems with bacteria growing prevented us to collect shedding data in the first 10 days post infection, consequently, we missed to identify the peak of shedding. Model simulations place the peak right at the start of the trial, however, this should be taken with caution, since there are no shedding data to train the model during the first 10 days (Figure 7). Therefore, while we still report the average trend of shedding during this period this has to be validated in the laboratory.

Figure 7

Download asset Open asset

Estimated dynamics of B. *bronchiseptica* shedding by time and type of infection.

The empirical shedding events (gray points), the estimated individual trajectories (smoothed thin black lines) and the estimated median group trends (blue lines), with the related 95%CIs (blue shadows), are reported. Level of shedding is presented as total daily event to scale up with model dynamics performed at one-day time step. For the B alone rabbits, we only report the group trend. Here, model prediction places the peak at the start of the trial (0 dpi), however, there were no shedding data to train the model during the first ten days and this result warns prudence.

Results from model fitting confirmed that bacteria shedding reached higher levels and lasted for a longer time in co-infected rabbits than in hosts where B. bronchiseptica was the only infection (Figure 7). For example, the peak value (CFU/day; 95% CIs) was: 14,185 (11.06–12,258,022) for BTG, 8,266 (5.89–10,119,754) for BG and 7,863 (29.57–2,213,310) for BT. For completeness we also report the peak value for B. bronchiseptica only rabbits, which was 678 CFU/day (168–2723); however, this is not indicative of the real value, as already noted above. Model simulations well captured the different trends of shedding among the co-infections. Rabbits from the BTG infection exhibited the most rapid increase of shedding, the highest levels and the slowest decline post-peak, suggesting a possible synergistic effect of the two helminths. The BT group showed the most rapid decline, while rabbits from the BG group maintained a trend intermediate between the BTG and BT groups (Figure 7).

Model fitting of individual time series described well the general trend and the high levels of shedding early in the infection (Figure 7 also compare it with Figure 2), particularly for BTG and BT. In contrast, the model fitted at the group level, showed a tendency to underestimate the average temporal profile, although we note that this is probably a visual effect. In fact, this trend reflected the large number of medium-low shedding events, including the fact that about 30% of these were classified as zero shedding and occurred throughout every co-infected group. The large variation in bacteria shed both within and between rabbits also contributed to this average pattern. Rabbits showed rapid changes in the level of shedding between consecutive samplings, where bursts of bacteria alternated to low or no shedding in a matter of a few days. This pattern was common among hosts, irrespective of their infection group. This rapid dynamic could not be explained by the intensity of infection and related immune responses at the scale used to record such data. Instead, it is possible that changes in local conditions, such as variation in the severity of tissue inflammation or control of bacteria at the mucosa level, including changes in the amount of mucus formation and expulsion, had a stronger impact on these fast changes. We also noted that co-infected rabbits sniffled more frequently than B. bronchiseptica only animals (Pathak’s pers. obs.) and, probably, further contributed to alter the frequency and intensity of these events.

Sensitivity analysis

To explore in more detail how variation in shedding between types of infection was related to changes in key immune parameters, a sensitivity analysis was performed on model performance (see Materials and methods and Appendix 2). Findings showed that an increase in the rate of bacterial neutralization, μ_c , and, secondly, in the mmune growth rate, μ_a, led to a proportional decrease of the peak of shedding and, to a lesser extent, the time to reach this peak (Appendix 2). Among the three immune variables, the strongest negative impact on shedding was caused by changes in neutrophil (c₁ or a₁) and IgG (c₃ or a₃) rates. As expected, an increase in bacterial growth rate, r, was associated with a higher peak of shedding and, consequently, a shorter time to reach this peak (Appendix 1—figures 1 and 2). Overall, the three co-infections exhibited comparable trends and were more sensitive to immune changes than B. bronchiseptica only infected rabbits. The weak response of the latter group was probably caused by the very low bacterial growth rate (r=0.06) and thus shedding, suggesting a stronger immune control.

Discussion

There is increasing evidence that gastrointestinal helminths can impact the severity and time course of respiratory bacterial infections (Brady et al., 1999; Lass et al., 2013; Diniz et al., 2010; Ezenwa et al., 2010; Babu and Nutman, 2016; Long et al., 2019); however, how they affect the dynamics of shedding remains to be determined. This is particularly important in regions where chronic helminthiasis co-circulate with respiratory pathogens that are endemic or cause seasonal outbreaks. This is also relevant in areas where antimicrobial resistance is emerging as a threat, such that focusing on the treatment of helminths could be a rapid and effective way to reduce onward bacterial transmission.

We investigated how two gastrointestinal-restricted helminth species, with contrasting infection dynamics, altered the shedding of the highly infectious B. bronchiseptica through the modulation of the immune response. Empirical findings showed that the two helminths, taken as a single species or as a pair, significantly increased the level, the variation and the duration of shedding, along with enhancing the frequency of supershedding events. Model estimates suggested that these changes were related to the significant growth of bacteria in the respiratory tract. The developed of an immune response controlled the infection and led to the waning but not the complete clearance of B. bronchiseptica. As a consequence of this, the level of shedding decreased with time but rabbits carried on shedding throughout the experiments. We found that changes in the relative contribution of neutrophils, IgA and IgG to bacterial neutralization could explain these dynamics and differences in the pattern of shedding between and within types of infection. However, the rapid changes in the magnitude of these events, commonly observed among rabbits, could not be explained by the relationship between infection and immunity at the scale of analysis selected for this study.

We modeled the dynamics of shedding as representative of the infection in the whole respiratory tract, yet, simulations provided a possible explanation of the contrasting trends in the lungs and nasal cavity. The decline of shedding in the first 90 days post infection could be explained by the immune-mediated clearance of B. bronchiseptica from the lungs, previously reported using data from these experiments (Pathak et al., 2010; Pathak et al., 2012; Thakar et al., 2012). In contrast, the reduced but ongoing shedding later in the experiments could be the result of the weak control of bacteria in the nasal cavity (Pathak et al., 2010; Pathak et al., 2012). We found that neutrophils quickly increased and contributed to bacteria neutralization early on, these were then followed by antibodies that developed more slowly and with an initial delay in the co-infected hosts. The two helminth species appeared not to change the general trend of these variables but the timing and intensity of their relative contribution and, with this, temporal changes in the level of shedding. We previously showed that both helminths elicit a Th-2 immune response (Murphy et al., 2011; Murphy et al., 2013; Cattadori et al., 2019), however, while T. retortaeformis was cleared, or significantly reduced from the small intestine, G. strigosum persisted with high intensities in the stomach (Cattadori et al., 2005; Cattadori et al., 2008; Cattadori et al., 2019; Murphy et al., 2011; Murphy et al., 2013; Mignatti et al., 2016). It is conceivable to think that the faster decline of shedding in rabbits from the BT than the BG or BTG group could have been caused by the stronger immune control of B. bronchiseptica, probably facilitated by the rapid clearance of T. retortaeformis (Thakar et al., 2012). The stronger immune control could also explain the tendency to the lower shedding and the fewer supershedding events observed. In contrast, the persistence of G. strigosum in the stomach probably contributed to delay the decline of shedding and to increase the magnitude of these events, including increasing the frequency of supershedding.

Interestingly, in the BTG triple infection the two helminths appeared to have a positive synergistic effect particularly in delaying the decline of shedding. Given these patterns, and considering that T. retortaeformis and G. strigosum are commonly found in populations of rabbits, it is not surprising if these helminths play an important role to the rapid transmission and persistence of B. bronchiseptica in natural settings. For example, the serological study of a free-living population of rabbits showed that the annual prevalence of B. bronchiseptica ranged between 88% and 97% (Pathak et al., 2010). Of these rabbits 65% were co-infected with B. bronchiseptica and G. strigosum, and both bacterial prevalence and helminth intensity increased with rabbit age. Although the study did not examine these patterns in dual or triple infections with T. retortaeformis, about 60% of adult rabbits were found to carry BTG infections (unpubl. data), suggesting that rabbits co-infected with helminths represent a large fraction of this host population and, together with the dual infections, are probably responsible for the high prevalence of B. bronchiseptica throughout the years.

Recent Boolean modeling of the whole immune network to B. bronchiseptica infection indicated that the combined effect of antigen-antibody complex and phagocytosis, activated by neutrophils among others, were important for bacterial removal from the lungs (Thakar et al., 2012). The model framework proposed in the current work adds to previous studies by linking a much simplified immune response to the dynamics of shedding through the intensity of infection in the whole respiratory tract. Model estimates showed that the rate of bacterial neutralization by IgA and IgG was stronger in rabbits from the BT group while IgA was particularly low in the BTG rabbits and, secondly, animals from the BG group. These differences could contribute to explain the fast reduction of shedding in the first group and the slower trend in the latter two groups. We found that the much lower rates of antibody neutralization in the B. bronchiseptica only rabbits were associated with very low bacterial growth and level of shedding, a pattern consistent with a strong immune control rather than a weak response. The important role of antibodies to B. bronchiseptica control was empirically highlighted in the mouse model where serum antibodies (which include IgG isotypes and IgA) cleared B. bronchiseptica from lungs and trachea of wild-type and B-cell-deficient mice within 3 days post inoculation (Kirimanjeswara et al., 2003). Similarly, IgA was found to be necessary for the reduction of bacterial numbers in the nasal cavity of mice (Wolfe et al., 2007).

Our model assumed that neutrophils were stimulated by, and directed against, B. bronchiseptica replication in the respiratory tract. The low bacterial shedding in the B group appeared to be associated with the prompt increase of neutrophils and their stronger rate of neutralization, compared to antibodies. In the co-infected rabbits, neutrophils showed to have a more variable role. The neutralization rate was consistently lower in rabbits from the BT group, appeared to have a mixed role with IgG in the BG group and dominated in the BTG group. In our previous work we found an early peak of neutrophils in the blood of rabbits infected with our two helminths and suggested that this was likely a consequence of gut bacteria infiltrating in the mucosa damaged during helminth establishment in the gastrointestinal tract (Murphy et al., 2011; Murphy et al., 2013). Although we were not able to disentangle the proportion of neutrophils that were stimulated by and directed against B. bronchiseptica, previous work by Thakar et al., 2012 found no temporal differences in bacterial clearance from the lungs of rabbits from the B and the BT groups, suggesting that there is probably no significant interference from the gut-stimulated neutrophils. The evidence that gut-restricted helminths could enhance the neutrophil response in the respiratory tract was showed in mice co-infected with Pseudomonas aeruginosa in the lungs and Heligmosomoides polygyrus in the gastrointestinal tract (Long et al., 2019). This study found a higher recruitment of neutrophils in the lungs together with an increase of CD4 +T cells and Th2 cytokine expression. Interestingly, survival rate was improved in co-infected mice than mice with bacteria alone, suggesting that neutrophils, rather than a Th2 response, dominated in the lungs. Work by Rolin et al., 2014 showed that TLR4-mutant mice infected with B. bronchiseptica poorly controlled bacterial growth and shedding. Animals exhibited heavy neutrophil infiltration in the lungs, however, the depletion of neutrophils did not affect the level of infection but decreased individual shedding. Taken together, we propose that the possible interference with a Th2 immune signal stimulated by T. retortaeformis and G. strigosum contributed to increase the level and duration of B. bronchiseptica shedding, although this anti-inflammatory response was not sufficiently strong to prevent bacteria clearance from the lower respiratory tract. We also note that while we focused on neutrophils and antibodies, this interference most likely involved other immune variables important in the network of immune reactions (Thakar et al., 2012).

A common feature among rabbits, irrespective of the infection group they belong to, was the rapid fluctuation in the level of bacteria shed over time, including the unpredicted supershedding events. We could not explain these rapid changes based on our scale of analysis. Most likely, other processes generated the exponential bursts and timing of (super)shedding observed. These could include, but are not limited to, immune-mediated control of bacteria turnover in the mucosa or stochastic events during bacteria release, such as random mucus discharge, and also variation in host behavior. Our finding is consistent with other systems where fast, intermittent shedding was reported both from single (Hadinoto et al., 2009; Schiffer et al., 2009; Rolin et al., 2014; Spencer et al., 2015; Slater et al., 2016) and co-infected hosts (Byrne et al., 2019; Kao et al., 2007). A compelling study by Hadinoto et al., 2009 on Epstein-Barr virus shedding in healthy carriers showed that the virus was rapidly and continuously shed in the saliva but the process of virus production and control was regulated by multiple factors, including the immune response. Importantly, they showed that variation in the level and frequency of shedding was a process that occurred at the individual level over time. Our study further supports this finding by showing variation in B. bronchiseptica shedding at three levels: between types of infection, between hosts with each infection type, and within each host over time. Our model framework, and the scale of our analysis, was able to describe the first two types of variation but prevented us from capturing the local mechanism responsible for the rapid changes observed at the individual level. More work at the local tissue level is needed to explain these patterns, including the processes generating supershedding events and the relative contribution of the lower and upper respiratory tract.

As illustrated by our findings, the definition of supershedder can become problematic when hosts are infected with more than one parasite/pathogen. Our analysis showed that a cut-off value based on single infection can lead to the overestimation of supershedding in co-infected rabbits. However, by pulling the data together we can calculate a common percentile threshold and provide a more accurate estimation. We used two closely related helminth species and if we consider the 99th percentile threshold, we found that between 18% and 32% of the co-infected rabbits generated at least one supershedding event, which greatly contrast with the lack of cases from the B. bronchiseptica only rabbits. These are significant percentages in our co-infected groups, especially if we consider that about 30% of the shedding events were null. Therefore, we can speculate that given the magnitude and frequency of supershedding, the risk of onwards transmission following a contact is likely to non-linearly increase for the co-infected hosts.

This study adds novel insights into the role of gastrointestinal helminths to the dynamics of B. bronchiseptica shedding. We used a parsimonious mechanism of immune regulation based on previous work by others and ourselves. While we focused on antibodies and neutrophils other immune variables could have contributed to the dynamics of shedding observed. Our model framework can be adapted to include the impact of other factors, such as components of the innate immune response or hierarchical relationships between immune variables. To reduce model complexity, we did not explicitly quantify the dynamics of infection of the two helminths although this can be explored in future work, including the role of B. bronchiseptica to the dynamics of shedding by the two helminths.

Given that one quarter of the global human population is infected with helminths and considering that respiratory infections are among the top 10 causes of death by infectious diseases worldwide, understanding the modulatory role of helminth species to respiratory infections is important for developing treatments targeted to specific co-infection settings. The ability to detect and, ideally, control the high shedders and/or supershedders is also critical for reducing the risk of disease outbreak and spread, and should not be overlooked any longer.

Materials and methods

Ethic statement

Request a detailed protocol

Animals were housed in individual cages with food and water ad libitum and a 12 hr day/night cycle, in compliance with Animal Welfare Act regulations as well as the Guide for the Care and Use of Laboratory Animals. All animal procedures, including infections with B. bronchiseptica and the two helminth species, weekly blood collection and pathogen/parasite sampling at fixed time, were approved by the Institutional Animal Care and Use Committee of The Pennsylvania State University (IACUC 26082). All animal work complied with guidelines as reported in the Guide for the Care and Use of Laboratory Animals, 8th ed. National Research Council of the National Academies, National Academies Press Washington DC.

Bacteria strain and culture

Request a detailed protocol

We used B. bronchiseptica strain RB50 for all experiments. Bacteria were grown on Bordet-Gengou (BG) agar supplemented with 10% defibrinated sheep blood and streptomycin (20 μg/ml). The inoculum was prepared by growing the bacteria in Stainer-Scholte (SS) liquid culture medium at 37 °C overnight. For the infection, bacteria were re-suspended in sterile phosphate-buffered saline (PBS) at a density of $5 * 10^{4}$ CFU/ml, which was confirmed by plating serial dilutions of the inoculum on BG blood agar plates in triplicate (Pathak et al., 2010).

Laboratory infections

Request a detailed protocol

B. bronchiseptica single infection (B) and dual infections with either T. retortaeformis (BT) or G. strigosum (BG) are described in detail in Pathak et al., 2010; Pathak et al., 2012 and Thakar et al., 2012. The triple infection (BTG) followed the same experimental design and laboratory procedures of the dual infections. Here, we describe the general design of the experiments. New Zealand White, two months old, male rabbits (Harlan, USA), were challenged with B. bronchiseptica and two helminth species as follow: i- B. bronchiseptica (B) single infection: 32 infected and 16 controls, ii- B. bronchiseptica-Graphidium strigosum (BG): 31 infected and 16 controls, iii- B. bronchiseptica-T. retortaeformis (BT): 32 infected and 16 controls, and iv- the three agents together (BTG): 32 infected and 16 controls. Infection was performed by pipetting in each nare $2.5 * 10^{4}$ CFU/ml of bacteria diluted in 0.5 ml of PBS. For the co-infections animals also received, simultaneously by gavage, a single inoculum of water (5 ml) with either 5500T. retortaeformis or 650 G. strigosum third stage infective larvae, or both. Helminth doses followed natural infections (Cattadori et al., 2005; Cattadori et al., 2008). Control animals were sham inoculated with 1 ml of sterile PBS in the nares and gavaged with 5 ml of water. The dynamics of infection and related immune responses were then followed for 120 days (150 days for B) by sacrificing four infected and two control animals at fixed days post infection, as follow, B: 3, 7, 15, 30, 59, 90, 120, 150 days; BG: 7, 14, 30, 44, 62, 76, 90, 120 days; BT: 5, 8, 15, 31, 47, 61, 91, 120 days; BTG 7, 14, 30, 45, 60, 75, 90, 120 days. This sampling design follows important time points in the life cycle of the three infections (Murphy et al., 2011; Murphy et al., 2013; Pathak et al., 2010), a shift of 1 or 2 days between infections was necessary a few times to adjust with our laboratory activities.

Bacteria shed enumeration

Request a detailed protocol

At the start of each infection, a subset of infected rabbits was randomly selected (i.e. B=16, BG = 23, BT = 20 and BTG = 24) to quantify the level of bacteria shed by contact with a BG-blood agar petri dish over time (Pathak et al., 2010). Shedding by contact with a surface mimics the natural transmission of a respiratory pathogen without disruption of the bacteria population through swabbing. Shedding was assessed once a week for B single infection and 2 (BT) or 3 (BG and BTG) times a week for the co-infected rabbits. The use of a different number of hosts and frequency of sampling was determined by logistical constraints (i.e. personnel availability). For the B. bronchiseptica group the use of a different substrate (plastic balls) the first 10 days post-infection was not successful and we missed to record the dynamics of shedding in those early days. At every sampling point, rabbits were allowed to interact with the petri dish by direct oral-nasal contact and for a maximum of 10 min. Plates were removed earlier if animals chewed the plastic or the agar, and the duration of each interaction was recorded. Plates were then incubated at 37 °C for 48 hr and colonies counted and scaled to the interaction time (CFU/s). If there was an interaction but plates resulted negative, shedding was considered to be null while the lack of interaction (e.g. animals were not interested in the plate) was recorded as a missing point.

Systemic immune response

Request a detailed protocol

As representative of the immune response to B. bronchiseptica, we selected neutrophils and species-specific antibodies IgA and IgG. Previous laboratory experiments and related modeling studies suggested that these three variables contribute to clear B. bronchiseptica from the lungs and to reduce the colonization in the nasal cavity (Orndorff et al., 1999; Kirimanjeswara et al., 2003; Wolfe et al., 2007; Thakar et al., 2007; Thakar et al., 2012). For example, results from modeling the immune network in the lungs of mice and rabbits showed that the lack of antibody production, by B cells deletion, prevented bacterial clearance (Thakar et al., 2007; Thakar et al., 2012). Similarly, peripheral neutrophils recruited via pro-inflammatory cytokines contributed to the activation of phagocytic cells and bacterial neutralization (Thakar et al., 2007; Thakar et al., 2012). Consistent with these findings, experimental studies showed that adoptive transfer of serum antibodies in mice led to the removal of B. bronchiseptica by day 3 post inoculation Kirimanjeswara et al., 2003. Methodologies to quantify neutrophils, IgA and IgG are described in Pathak et al., 2010; Pathak et al., 2012 and Thakar et al., 2012. Briefly, for every rabbit blood was collected once a week for neutrophils and twice a week for antibodies. Neutrophil concentration was measured using whole blood (0.2 ml) stored in EDTA (Sartorius, Germany) and analyzed using Hemavet-3 hematology system (Drew Scientific, USA). Species-specific IgA and IgG were estimated from blood serum using B. bronchiseptica as a source of antigen and ELISA (Pathak et al., 2010). Measurements were performed in duplicates with all plates having high, low and background controls. Values were expressed as immunosorbent Optical Densities (OD) and then standardized into Optical Density index (Murphy et al., 2011). Plate preparation and dilutions, including the preparation of high (strongly reacting animals) and low (non-reacting animals from prior to the infection) pools and checkerboard titrations, are detailed elsewhere (Pathak et al., 2010; Pathak et al., 2012; Thakar et al., 2012). For the triple infection, antibody quantification and dilutions followed Pathak et al., 2010. The baseline profile of the three immune variables were available from blood data collected from every rabbit the week before the infection; if individual data were missing we used the average value from control rabbits. Neutrophils rapidly returned to baseline levels following the initial infection and the equilibrium value was available as the average from control rabbits sampled throughout the experiments, specifically cells/ml B: 1.323, BG: 2.730, BT: 1.819 and BTG: 2.289.

Model framework

Request a detailed protocol

The within-host mechanisms that affect the dynamics of B. bronchiseptica shedding were examined in two steps. First, we developed a dynamical model that describes the dynamics of neutrophils, specific IgA and IgG, and their interaction with the bacterial infection in the whole respiratory tract. Second, a Bayesian approach was then used to link this dynamical model to the empirical longitudinal data by: i- fitting the model to the three immune variables of every rabbit and ii- fitting the intensity of B. bronchiseptica infection, estimated from the dynamical model, to the experimental shedding data. Below we describe the dynamical model and in the next sections we report on the observation model and model fitting to the empirical data, including how parameter calibration was performed. For completeness, we also include two additional sections: i- model validation that explores the accuracy of our model and ii-model selection where we compare three different model formulations.

Dynamical model

Request a detailed protocol

Different modeling approaches have been applied to study the within-host dynamics of infection of Bordetella species and related immune response. Work by de Graaf et al., 2014 on B. pertussis followed a phenomenological approach based on a Bayesian hierarchical framework that described the rise and decline of IgG data as a response to the infection. In contrast, previous studies by colleagues and us focused on a network-based dynamical model wherein relationships between components of the immune system against B. bronchiseptica or B. pertussis infection were examined using experimental data and Boolean transfer functions (Thakar et al., 2007; Thakar et al., 2012). This approach allowed us to follow the complexity of the immune system both at the systemic and localized level of infection and the consequences on bacterial control. In the current study, we developed a deterministic dynamical model that explored the relationship between bacterial shedding and host immunity using an individual based Bayesian approach and laboratory data. Here, we simplified the immune response down to three variables and examined their direct effect on the dynamics of infection and related shedding. The time-dependent interactions between immune variables and B. bronchiseptica infection are described by the following system of ordinary differential equations:

\frac{d x_{1, i}}{d t} = a_{1, i} y_{i} - b_{1, i} (x_{1, i} - x_{1}^{*})

\frac{d x_{2, i}}{d t} = a_{2, i} y_{i} - b_{2} x_{2, i}

\frac{d x_{3, i}}{d t} = a_{3, i} y_{i} - b_{3} x_{3, i}

\frac{d y_{i}}{d t} = (r_{i} - c_{1, i} x_{1, i} - c_{2} x_{2, i} - c_{3} x_{3, i}) * y_{i}

where $x_{1, i} (t)$ , $x_{2, i} (t)$ and $x_{3, i} (t)$ are neutrophils, IgA, and IgG, respectively, of rabbit i^th at time t. $y_{i} (t)$ is the intensity of bacteria in the respiratory tract of this same rabbit at time t. The parameters $a_{1, i}$ , $b_{1, i}$ and $c_{1, i}$ describe the rate of growth, decay and bacteria clearance, respectively, for neutrophils in rabbit i^th. Similarly, $a_{2, i}, b_{2}$ and $c_{2, i}$ represent the per capita rates of growth, decay and bacterial clearance for IgA, while $a_{3, i}$ , b₃ and $c_{3, i}$ are the rates representing IgG (Table 1). The baseline immune conditions of the host before the infection are $x_{1, i} (0)$ , $x_{2, i} (0)$ and $x_{3, i} (0)$ , while $x_{1}^{*}$ describes the equilibrium level of neutrophils post infection; this was not included for IgA or IgG since the time to reach equilibrium extended beyond the course of the experiments. The parameter r_i represents the per capita bacterial replication rate. The full description of model parameters and variables is reported in Table 1. In our model formulation, we made the parsimonious assumption that the activation and response of neutrophils and specific IgA and IgG proportionately increase with the intensity of infection, $y_{i} (t)$ , in the whole respiratory tract. Similarly, bacterial neutralization occurs through the additive effect of these three immune variables. These direct interactions are a large simplification of a more complex immune process (Thakar et al., 2007; Thakar et al., 2012), and follow a classical Lotka-Volterra type of relationship commonly used to describe within-host processes of infection (Mohtashemi and Levins, 2001; Pugliese and Gandolfi, 2008; Fenton and Perkins, 2010; de Graaf et al., 2014; Vanalli et al., 2020). To reduce model complexity we assumed that the two helminths affect B. bronchiseptica infection, and thus shedding, by altering the magnitude and time course of neutrophil, IgA and IgG responses.

Observation model

Request a detailed protocol

We applied the dynamical model to every rabbit (except for the B. bronchiseptica group). Individual parameters could vary among hosts as independent samples from a joint log-normal distribution, while parameters that were shared among rabbits (see below) were kept the same for that group. This hierarchical set-up allowed us to have individual responses that varied in term of amplitude, time to peak and decay rate, while keeping these trends from deviating too much from each other. For a given rabbit i^th at time t, every empirical immune variable, log-transformed, was assumed to follow a normal distribution with mean, µ, and variance, $σ$ , as:

\log (N e_{i, t}) \sim N (\log (x_{1, i} (t)), σ_{x_{1}})

\log (I_{A, i, t}) \sim N (\log (x_{2, i} (t)), σ_{x_{2}})

\log (I_{G, i, t}) \sim N (\log (x_{3, i} (t)), σ_{x_{3}})

The empirical amount of B. bronchiseptica shed by rabbit i^th at time t, is directly proportional to the level of infection $y_{i} (t)$ and is assumed to be representative of the intensity of infection in the whole respiratory tract. The probability of having a shedding event is independent of time since inoculation, in that shedding can occur anytime during the experiment and anytime during the interaction with the petri dish. Shedding was then related to the dynamics of infection via a zero-inflated log-normal relationship, to account for the high fraction of events in which rabbits did not shed (i.e. null shedding) despite been infected (Figure 3), as:

\log (S_{i, t}) \sim w_{0, i} + (1 - w_{0, i}) * N (\log (y_{i} (t)), σ_{y})

where $w_{0, i}$ is the fraction of events when a given host did not shed. For consistency with the immune data, collected on a daily scale, the estimated shedding was then scaled up and quantified as total amount of bacteria shed by a rabbit in a day.

Parameter calibration and model fitting

Request a detailed protocol

The dynamical model described in equations (1)-(4) includes parameters that are shared among rabbits within the same group, specifically:

θ_{1} = {μ_{a_{1}}, σ_{a_{1}}, μ_{b_{1}}, σ_{b_{1}}, μ_{c_{1}}, σ_{c_{1}}, μ_{a_{2}}, σ_{a_{2}}, b_{2}, μ_{c_{2}}, σ_{c_{2}}, μ_{a_{3}}, σ_{a_{3}}, b_{3}, μ_{c_{3}}, σ_{c_{3}}, μ_{r}, σ_{r}

where $μ$ and $σ$ are the means and standard deviations, respectively, of the corresponding normal distributions of the parameters $a$ , $b$ , $c$ and $r$ for each group. Likewise, the set of parameters $θ_{2, i}$ represents the estimates of every individual rabbit within a group, as:

θ_{2, i} = {a_{1, i}, b_{1, i}, c_{1, i}, a_{2, i}, c_{2, i}, a_{3, i}, c_{3, i}, r_{i}, w_{0, i}}_{(i \in 1, . . ., N)}

where N is the number of infected rabbits in each group. The resulting hyper prior distribution for rabbit i^th is:

p (θ_{2, i}) = ϕ (l o g (a_{1, i}) | μ_{a_{1}}, σ_{a_{1}}) * ϕ (l o g (b_{1, i}) | μ_{b_{1}}, σ_{b_{1}}) * ϕ (l o g (c_{1, i}) | μ_{c_{1}}, σ_{c_{1}}) * ϕ (l o g (a_{2, i}) | μ_{a_{2}}, σ_{a_{2}}) * ϕ (l o g (c_{2, i}) | μ_{c_{2}}, σ_{c_{2}}) * ϕ (l o g (a_{3, i}) | μ_{a_{3}}, σ_{a_{3}}) * ϕ (l o g (c_{3, i}) | μ_{c_{3}}, σ_{c_{3}}) * ϕ (l o g (r_{i}) | μ_{r}, σ_{r})

The resulting joint-likelihood function for the observed longitudinal time-series of all infected rabbits is:

L (θ_{1}, θ_{2}) \propto p (θ_{1}) * \prod_{i = 1}^{N} p (θ_{2, i}) * \prod_{i = 1}^{N} \prod_{j = 1}^{n_{i}} ϕ (l o g (N e_{i, j}) | l o g (x_{1, i, j}), σ_{x_{1}}) * \prod_{i = 1}^{N} \prod_{k = 1}^{n_{1, i}} ϕ (l o g (I_{A, i, k}) | l o g (x_{2, i, k}), σ_{x_{2}}) * \prod_{i = 1}^{N} \prod_{l = 1}^{n_{2, i}} ϕ (l o g (I_{G, i, l}) | l o g (x_{3, i, l}), σ_{x_{3}}) * \prod_{i = 1}^{N} \prod_{h = 1}^{n_{S, i}} ϕ (l o g (S_{i, h}) | l o g (y_{i, h}), σ_{y})

where n_i, $n_{1, i}$ , $n_{2, i}$ , and $n_{S, i}$ are the number of measurements for neutrophils, IgA, IgG and bacteria shed, respectively, for rabbit i^th. The four subscripts j, k, l, h are the time index for neutrophils, IgA, IgG and bacteria shed, respectively, at the specific time points for rabbit i^th. The notation $ϕ (x | μ, σ^{2})$ represents the normal probability density of $x$ with mean µ and variance $σ^{2}$ . Since there is no knowledge on the values of $θ_{1}$ we used normal priors as reported in Table 2.

Model fitting, and the generation of posterior distributions for the parameters estimated, was performed using Hamiltonian Monte Carlo (HMC) algorithm implemented in Stan package version 2.18.0. Briefly, the algorithm was implemented as follow:

1. Start sampling from the joint normal prior distributions for $θ$ , here referring to all parameters included in $θ_{1}$ and $θ_{2}$ , (the initial values $θ^{0}$ could be user-specific or specified by Stan). Set the initial value of $θ$ as the current $θ^{*}$ .

2. Evaluate momentum vector for the current $θ^{*}$ (Note that in HMC, the momentum vector $ℒ (θ)$ is described by its ‘kinetic’ and ‘potential’ energies).

3. The HMC proposes a new $θ^{*}$ by sampling from the posterior distributions of $p (θ | y)$ with a given stepsize which is optimized by the algorithm.

4. To account for numerical errors during integration, a Metropolis acceptance step is applied. If the HMC proposed momentum vector has a higher probability than the previous parameter values, the proposed values will be updated to the current values and used to initialize the next iteration. If the Metropolis acceptant rejects the proposed values, the previous values are returned and used to initialize the next iteration. Repeat the following steps 2–4 until the maximum number of iterations is reached. This procedure is for a single chain, there are four parallel chains, each with 100,000 iterations including 40,000 burn-in iterations. This procedure was repeated for every rabbit in each of the four types of infection.

5. The generated sampling distribution represents the posterior distribution for each parameter, once the warm-up iterations have been removed. The posterior distribution is now used to quantify the statistics of interest (e.g. mean, median, 95% CI) for each group.

Model simulations were performed at one-day time step increment and parameters were scaled to account for differences in their relative magnitude. To provide meaningful results and to avoid problems with model convergence, rabbit time series with three, or less, shedding points larger than zero were excluded. This led us to used: 17 BT, 22 BG, 23 BTG, and 7 B rabbits. Given the small number of rabbits, and individual time series with few shedding events larger than zero, the modeling of the B. bronchiseptica only rabbits was carried out only at the group level.

To ensure proper model convergence, we assessed: i- the scale reduction factor $\hat{R}$ on split chains, to confirm that the value is close to 1 for each parameter, a $\hat{R} > 1.1$ is usually an indicator of a fit problem, and ii- the crude measure of the effective sample size $n_{e f f}$ ,this value should be close to 1 for a good parameter estimate. Moreover, the Rstan package allows us to set the maximum trajectory length to avoid infinite loops that can occur for non-identified models. If a high proportion of the interactions saturate the maximum threshold, the model is not effectively sampled. Finally, we checked the number of divergences, the trace plot of the Markov Chain time-series and the plot of the priors and posteriors for each parameter. All parameter estimates were diagnosed with the above criteria to ensure that the chains were well-mixed with no divergent transitions post-burning (Appendix 3).

Sensitivity analysis

Request a detailed protocol

To examine how key immune parameters affected shedding, a sensitivity analysis was performed (Appendix 2). As representative of shedding, we considered the peak of shedding and the time to reach the peak, including bacteria growth rate r, while for the immune variables we selected the rate of bacterial neutralization $μ_{c}$ and the growth rate $μ_{a}$ . Analyses were carried out using mean values for each type of infection. Specifically, for either c or a the estimated optimal value presented in Table 2 was changed by an incremental percentage, while keeping all the other parameters fixed, and the consequent changes in the amplitude and timing of the shedding peak was quantified. The same exercise was repeated using incremental changes of r; 1000 simulations were performed for each parameter at each step.

Model validation

Request a detailed protocol

As part of our technical check of model performance, to investigate whether the Bayesian model could recover the starting parameter values we performed a simulation-based validation using parameters at the group level. Briefly, we randomly drew $n$ sets of parameters $θ$ from our prior distributions (n=1,000 iterations) and generated a y_i for each parameter set of $θ_{i}$ from our dynamical model. We then fitted the simulated dataset to our Bayesian framework to confirm that the posterior distributions included the known sets of parameters. The model showed a good ability to recover the original parameters (Appendix 4).

Model selection

Request a detailed protocol

To examine the effect of different immune responses on the pattern of shedding, three model formulations were compared: i- neutrophils, IgA and IgG (full model), ii: neutrophils and specific IgA (reduced model A) and iii: neutrophils and specific IgG (reduced model G), where neutrophils were kept as the fundamental variable. To simplify the workflow, model selection was performed using data from each group in the form of time series of geometric means for every variable for each of the four infections, group level parameters were also used. The best model was selected based on the best compromise between goodness of fit and parsimony, according to the Bayesian information criterion (BIC). Results showed that the full model was statistically better in capturing the dynamics of B. bronchiseptica shedding, than the reduced formulations (Appendix 5). This model was then used to perform simulations both at the level of individual rabbit and infection group, with the exception of the B. bronchiseptica only rabbits that were examined at the group level.

Appendix 1

Parameter estimates from posterior distributions

Averages of posterior parameter estimations at the individual and group level were examined for the three co-infection groups (BG, BT and BTG); B. bronchiseptica alone rabbits were excluded since simulations were based on the whole group. The decay rates of IgA and IgG, b₂ and b₃, respectively, were assumed to be the same for rabbits of each group and were not reported. Results show clear rabbit variation in B. bronchiseptica dynamics of infection and immune responses, both within and between groups (Appendix 1—figures 1–3).

Appendix 1—figure 1

Download asset Open asset

Parameter estimates from posterior distribution for BG.

Individual average (black points) with 95% CIs (black segments), and group average (named ’Pop’, red points) with 95% CIs (red segments), are reported. X-axis lists the rabbit’s ID every 5 animals using the same left-to-right ID order, and host sampling scheme, as detailed in Figure 1 of the main text.

Appendix 1—figure 2

Download asset Open asset

Parameter estimates from posterior distribution for BT.

Full details in Appendix 1—figure 1.

Appendix 1—figure 3

Download asset Open asset

Parameter estimates from posterior distribution for BTG.

Full details in Appendix 1—figure 1.

Appendix 2

Sensitivity analysis

Sensitivity analysis showed that the three co-infections exhibited comparable trends and were more sensitive to immune changes than B. bronchiseptica only rabbits (Appendix 2—figures 1–2); additional details in the main text.

Appendix 2—figure 1

Download asset Open asset

Relationship between B.bronchiseptica peak of shedding and percentile changes in neutrophils, IgA, and IgG.

The growth a, and neutralizing c, rates of neutrophils (subscript 1), IgA (subscript 2) and IgG (subscript 3), including bacterial growth rate r, are reported by group. Mean estimates with 95% CIs are presented.

Appendix 2—figure 2

Download asset Open asset

Relationship between B.bronchiseptica time to peak of shedding and percentile changes in neutrophils, IgA, and IgG.

Full details in Appendix 2—figure 1.

Appendix 3

Model convergence

Here, we present convergence plots for the four model chains for $θ_{1 - 10}$ , where each chain was run with 100,000 iterations (Appendix 3—figures 1–4). For each θ and infection group the chains were well mixed and largely overlapping. For completeness, we have also included the correlations and related plots from the posterior distributions (Appendix 3—figures 5–8) and the scale reduction factors for each parameter (Appendix 3—table 1). Additional details in Materials and methods.

Appendix 3—figure 1

Download asset Open asset

MCMC trace plot for B.

Different colors represent the four chains.

Appendix 3—figure 2

Download asset Open asset

MCMC trace plot for BG.

Different colors represent the four chains.

Appendix 3—figure 3

Download asset Open asset

MCMC trace plot for BT.

Different colors represent the four chains.

Appendix 3—figure 4

Download asset Open asset

MCMC trace plot for BTG.

Different colors represent the four chains.

Appendix 3—figure 5

Download asset Open asset

MCMC posteriors from B.

Histograms (diagonal), pair-wise scatterplots (upper right) and related correlations (lower left) from the parameter posteriors.

Appendix 3—figure 6

Download asset Open asset

MCMC posteriors from BG.

Histograms (diagonal), pair-wise scatterplots (upper right) and related correlations (lower left) from the parameter posteriors.

Appendix 3—figure 7

Download asset Open asset

MCMC posteriors from BT.

Histograms (diagonal), pair-wise scatterplots (upper right) and related correlations (lower left) from the parameter posteriors.

Appendix 3—figure 8

Download asset Open asset

MCMC posteriors from BTG.

Histograms (diagonal), pair-wise scatterplots (upper right) and related correlations (lower left) from the parameter posteriors.

Appendix 3—table 1

Scale reduction factors for each parameter of the four types of infection.

Parameters	B	BG	BT	BTG
$μ_{a 1}$	1.000234	1.000103	1.000252	1.000047
$μ_{b 1}$	1.000166	1.000205	1.000122	1.000073
$μ_{c 1}$	1.000053	1.000047	0.999977	1.000423
$μ_{a 2}$	1.000120	1.000174	1.000285	0.999986
$μ_{b 2}$	1.000043	1.000049	1.000051	1.000128
$μ_{c 2}$	1.000060	0.999973	1.000240	1.000316
$μ_{a 3}$	1.000068	1.000166	1.000027	1.000059
$μ_{b 3}$	1.000013	0.999927	0.999919	1.000089
$μ_{c 3}$	0.999967	1.000016	1.000224	1.000204
$μ_{r}$	0.999972	0.999961	1.000139	1.000517
log(LLH)	1.000650	1.000406	1.000521	1.000401

Appendix 4

Model validation

We investigated the ability of the model to recover the correct parameter values and Appendix 4—figure 1 shows good model performance and accuracy for every estimated parameter. Additional details in Materials and methods.

Appendix 4—figure 1

Download asset Open asset

Model validation using group-level parameters.

The vertical lines represent the true parameters that were used to simulate the datasets.

Appendix 5

Model selection

Three alternative dynamic models were tested: i- full model (neutrophils +IgA + IgG), ii- reduced model A (neutrophils +IgA) and iii: reduce model G (neutrophils +IgG); neutrophils were kept as the fundamental variable (Appendix 5—table 1). Additional details in Materials and methods.

Appendix 5—table 1

Bayesian Information Criterion (BIC) and level of model complexity (N), including pair-wise Δ BIC with the best fitted model.

Models	BIC-ΔBIC	B	BG	BT	BTG	N
Neutrophils+IgA+IgG	BIC	-652.9	-351.2	594.8	-610.5	12
Neutrophils+IgA	BIC	-251.3	-302.2	-570.9	-507.7	9
Neutrophils+IgA	ΔBIC	401.5	48.9	23.9	102.8
Neutrophils+IgG	BIC	-314.8	-105.6	-442.0	-423.9	9
Neutrophils+IgG	ΔBIC	338.1	245.6	152.7	186.5

Data availability

All empirical data used in this study are available at https://doi.org/10.5061/dryad.g79cnp5sx. File titles: Nguyen-Pathak-Cattadori eLife Bordetella shedding.csv; Nguyen-Pathak-Cattadori eLife neutrophils.csv; Nguyen-Pathak-Cattadori eLife IgA ODI.csv; Nguyen-Pathak-Cattadori eLife IgG ODI.csv; README.doc.

The following data sets were generated

(2022) Dryad Digital Repository
Gastrointestinal Helminths Increase Bordetella bronchiseptica Shedding and Host Variation in Supershedding.
https://doi.org/10.5061/dryad.g79cnp5sx

References

1. Babu S
2. Nutman TB
(2016) Helminth-tuberculosis co-infection: an immunologic perspective
Trends in Immunology 37:597–607.
https://doi.org/10.1016/j.it.2016.07.005
- PubMed
- Google Scholar
(1999) Fasciola hepatica suppresses a protective Th1 response against Bordetella pertussis
Infection and Immunity 67:5372–5378.
https://doi.org/10.1128/IAI.67.10.5372-5378.1999
- PubMed
- Google Scholar
(2008) Coinfection of pigs with porcine respiratory coronavirus and Bordetella bronchiseptica
Veterinary Microbiology 128:36–47.
https://doi.org/10.1016/j.vetmic.2007.09.025
- PubMed
- Google Scholar
1. Byrne CM
2. Johnston C
3. Orem J
4. Okuku F
5. Huang ML
6. Selke S
7. Wald A
8. Corey L
9. Schiffer JT
10. Casper C
11. Coombs D
12. Gantt S
(2019) Increased oral epstein-barr virus shedding with HIV-1 co-infection is due to a combination of B cell activation and impaired cellular immune control
Pathology 11:587063.
https://doi.org/10.1101/587063
- Google Scholar
(2005) Peak shift and epidemiology in a seasonal host-nematode system
Proceedings. Biological Sciences 272:1163–1169.
https://doi.org/10.1098/rspb.2004.3050
- PubMed
- Google Scholar
(2008) Parasite co-infection and interaction as drivers of host heterogeneity
International Journal for Parasitology 38:371–380.
https://doi.org/10.1016/j.ijpara.2007.08.004
- PubMed
- Google Scholar
(2019) External disturbances impact helminth-host interactions by affecting dynamics of infection, parasite traits, and host immune responses
Ecology and Evolution 9:13495–13505.
https://doi.org/10.1002/ece3.5805
- PubMed
- Google Scholar
(2008) Super-shedding and the link between human infection and livestock carriage of Escherichia coli O157
Nature Reviews. Microbiology 6:904–912.
https://doi.org/10.1038/nrmicro2029
- PubMed
- Google Scholar
1. Chen WJ
2. Yang JY
3. Lin JH
4. Fann CSJ
5. Osyetrov V
6. King CC
7. Chen YMA
8. Chang HL
9. Kuo HW
10. Liao F
11. Ho MS
(2006) Nasopharyngeal shedding of severe acute respiratory syndrome-associated coronavirus is associated with genetic polymorphisms
Clinical Infectious Diseases 42:1561–1569.
https://doi.org/10.1086/503843
- PubMed
- Google Scholar
(2014) A two-phase within-host model for immune response and its application to serological profiles of pertussis
Epidemics 9:1–7.
https://doi.org/10.1016/j.epidem.2014.08.002
- PubMed
- Google Scholar
(2010) Presence of intestinal helminths decreases T helper type 1 responses in tuberculoid leprosy patients and may increase the risk for multi-bacillary leprosy
Clinical and Experimental Immunology 161:142–150.
https://doi.org/10.1111/j.1365-2249.2010.04164.x
- PubMed
- Google Scholar
(2018) The impact of past vaccination coverage and immunity on pertussis resurgence
Science Translational Medicine 10:eaaj1748.
https://doi.org/10.1126/scitranslmed.aaj1748
- PubMed
- Google Scholar
(2010) Hidden consequences of living in a wormy world: nematode‐induced immune suppression facilitates tuberculosis invasion in African buffalo
The American Naturalist 176:613–624.
https://doi.org/10.1086/656496
- PubMed
- Google Scholar
1. Fenton A
2. Perkins SE
(2010) Applying predator-prey theory to modelling immune-mediated, within-host interspecific parasite interactions
Parasitology 137:1027–1038.
https://doi.org/10.1017/S0031182009991788
- PubMed
- Google Scholar
1. Goodnow RA
(1980) Biology of Bordetella bronchiseptica
Microbiological Reviews 44:722–738.
https://doi.org/10.1128/mr.44.4.722-738.1980
- PubMed
- Google Scholar
(2012) Shedding light on Salmonella carriers
Trends in Microbiology 20:320–327.
https://doi.org/10.1016/j.tim.2012.04.004
- PubMed
- Google Scholar
1. Gopinath S
2. Hotson A
3. Johns J
4. Nolan G
5. Monack D
(2013) The systemic immune state of super-shedder mice is characterized by a unique neutrophil-dependent blunting of th1 responses
PLOS Pathogens 9:e1003408.
https://doi.org/10.1371/journal.ppat.1003408
- PubMed
- Google Scholar
(2007) Transmission consequences of coinfection: cytokines writ large?
Trends in Parasitology 23:284–291.
https://doi.org/10.1016/j.pt.2007.04.005
- PubMed
- Google Scholar
(2009) The dynamics of EBV shedding implicate a central role for epithelial cells in amplifying viral output
PLOS Pathogens 5:e1000496.
https://doi.org/10.1371/journal.ppat.1000496
- PubMed
- Google Scholar
(2018) Bordetella bronchiseptica colonization limits efficacy, but not immunogenicity, of live-attenuated influenza virus vaccine and enhances pathogenesis after influenza challenge
Frontiers in Immunology 9:2255.
https://doi.org/10.3389/fimmu.2018.02255
- PubMed
- Google Scholar
1. Kao RR
2. Gravenor MB
3. Charleston B
4. Hope JC
5. Martin M
6. Howard CJ
(2007) Mycobacterium bovis shedding patterns from experimentally infected calves and the effect of concurrent infection with bovine viral diarrhoea virus
Journal of the Royal Society, Interface 4:545–551.
https://doi.org/10.1098/rsif.2006.0190
- PubMed
- Google Scholar
(2001) Dynamics of the 2001 UK foot and mouth epidemic: stochastic dispersal in a heterogeneous landscape
Science 294:813–817.
https://doi.org/10.1126/science.1065973
- PubMed
- Google Scholar
(2003) Role of antibodies in immunity to Bordetella infections
Infection and Immunity 71:1719–1724.
https://doi.org/10.1128/IAI.71.4.1719-1724.2003
- PubMed
- Google Scholar
(2016) Association between pneumocystis spp. and co-infections with Bordetella bronchiseptica, mycoplasma hyopneumoniae and pasteurella multocida in austrian pigs with pneumonia
Veterinary Journal 207:177–179.
https://doi.org/10.1016/j.tvjl.2015.11.003
- PubMed
- Google Scholar
1. Lass S
2. Hudson PJ
3. Thakar J
4. Saric J
5. Harvill E
6. Albert R
7. Perkins SE
(2013) Generating super-shedders: co-infection increases bacterial load and egg production of a gastrointestinal helminth
Journal of the Royal Society, Interface 10:20120588.
https://doi.org/10.1098/rsif.2012.0588
- PubMed
- Google Scholar
1. Leung NHL
2. Xu C
3. Ip DKM
4. Cowling BJ
(2015) Review article: the fraction of influenza virus infections that are asymptomatic: a systematic review and meta-analysis
Epidemiology 26:862–872.
https://doi.org/10.1097/EDE.0000000000000340
- PubMed
- Google Scholar
(2005) Superspreading and the effect of individual variation on disease emergence
Nature 438:355–359.
https://doi.org/10.1038/nature04153
- PubMed
- Google Scholar
1. Long GH
2. Sinha D
3. Read AF
4. Pritt S
5. Kline B
6. Harvill ET
7. Hudson PJ
8. Bjørnstad ON
(2010) Identifying the age cohort responsible for transmission in a natural outbreak of Bordetella bronchiseptica
PLOS Pathogens 6:e1001224.
https://doi.org/10.1371/journal.ppat.1001224
- PubMed
- Google Scholar
1. Long SR
2. Lanter BB
3. Pazos MA
4. Mou H
5. Barrios J
6. Su CW
7. Wang ZQ
8. Walker WA
9. Hurley BP
10. Shi HN
(2019) Intestinal helminth infection enhances bacteria-induced recruitment of neutrophils to the airspace
Scientific Reports 9:1–13.
https://doi.org/10.1038/s41598-019-51991-3
- PubMed
- Google Scholar
(2010) Influenza virus coinfection with Bordetella bronchiseptica enhances bacterial colonization and host responses exacerbating pulmonary lesions
Microbial Pathogenesis 49:237–245.
https://doi.org/10.1016/j.micpath.2010.06.004
- PubMed
- Google Scholar
(2016) Host immunity shapes the impact of climate changes on the dynamics of parasite infections
PNAS 113:2970–2975.
https://doi.org/10.1073/pnas.1501193113
- PubMed
- Google Scholar
1. Mohtashemi M
2. Levins R
(2001) Transient dynamics and early diagnostics in infectious disease
Journal of Mathematical Biology 43:446–470.
https://doi.org/10.1007/s002850100103
- PubMed
- Google Scholar
(2011) Explaining patterns of infection in free-living populations using laboratory immune experiments
Parasite Immunology 33:287–302.
https://doi.org/10.1111/j.1365-3024.2011.01281.x
- PubMed
- Google Scholar
(2013) A co-infection with two gastrointestinal nematodes alters host immune responses and only partially parasite dynamics
Parasite Immunology 35:421–432.
https://doi.org/10.1111/pim.12045
- PubMed
- Google Scholar
1. Orndorff PE
2. Harvill ET
3. Cotter PA
4. Yuk MH
5. Miller JF
(1999) Probing the function of Bordetella bronchiseptica adenylate cyclase toxin by manipulating host immunity
Infection and Immunity 67:1493–1500.
https://doi.org/10.1128/IAI.67.3.1493-1500.1999
- Google Scholar
(2010) Immune regulation of a chronic bacteria infection and consequences for pathogen transmission
BMC Microbiology 10:1–9.
https://doi.org/10.1186/1471-2180-10-226
- PubMed
- Google Scholar
(2012) Immuno-epidemiology of chronic bacterial and helminth co-infections: observations from the field and evidence from the laboratory
International Journal for Parasitology 42:647–655.
https://doi.org/10.1016/j.ijpara.2012.04.011
- PubMed
- Google Scholar
1. Pugliese A
2. Gandolfi A
(2008) A simple model of pathogen-immune dynamics including specific and non-specific immunity
Mathematical Biosciences 214:73–80.
https://doi.org/10.1016/j.mbs.2008.04.004
- PubMed
- Google Scholar
(2014) TLR2 signaling decreases transmission of streptococcus pneumoniae by limiting bacterial shedding in an infant mouse influenza A co-infection model
PLOS Pathogens 10:e1004339.
https://doi.org/10.1371/journal.ppat.1004339
- PubMed
- Google Scholar
1. Rolin O
2. Smallridge W
3. Henry M
4. Goodfield L
5. Place D
6. Harvill ET
(2014) Toll-like receptor 4 limits transmission of Bordetella bronchiseptica
PLOS ONE 9:e85229.
https://doi.org/10.1371/journal.pone.0085229
- PubMed
- Google Scholar
1. Schiffer JT
2. Abu-Raddad L
3. Mark KE
4. Zhu J
5. Selke S
6. Magaret A
7. Wald A
8. Corey L
(2009) Frequent release of low amounts of herpes simplex virus from neurons: results of a mathematical model
Science Translational Medicine 1:7ra16.
https://doi.org/10.1126/scitranslmed.3000193
- PubMed
- Google Scholar
1. Schulz BS
2. Kurz S
3. Weber K
4. Balzer HJ
5. Hartmann K
(2014) Detection of respiratory viruses and Bordetella bronchiseptica in dogs with acute respiratory tract infections
Veterinary Journal 201:365–369.
https://doi.org/10.1016/j.tvjl.2014.04.019
- PubMed
- Google Scholar
1. Sheth PM
2. Danesh A
3. Sheung A
4. Rebbapragada A
5. Shahabi K
6. Kovacs C
7. Halpenny R
8. Tilley D
9. Mazzulli T
10. MacDonald K
11. Kelvin D
12. Kaul R
(2006) Disproportionately high semen shedding of HIV is associated with compartmentalized cytomegalovirus reactivation
The Journal of Infectious Diseases 193:45–48.
https://doi.org/10.1086/498576
- PubMed
- Google Scholar
1. Slater N
2. Mitchell RM
3. Whitlock RH
4. Fyock T
5. Pradhan AK
6. Knupfer E
7. Schukken YH
8. Louzoun Y
(2016) Impact of the shedding level on transmission of persistent infections in mycobacterium avium subspecies paratuberculosis (MAP)
Veterinary Research 47:1–12.
https://doi.org/10.1186/s13567-016-0323-3
- PubMed
- Google Scholar
(2013) Kinetics of coinfection with influenza A virus and streptococcus pneumoniae
PLOS Pathogens 9:e1003238.
https://doi.org/10.1371/journal.ppat.1003238
- PubMed
- Google Scholar
(2015) “ Super” or just “ above average”? supershedders and the transmission of Escherichia coli O157: H7 among feedlot cattle
Journal of the Royal Society, Interface 12:0446.
https://doi.org/10.1098/rsif.2015.0446
- PubMed
- Google Scholar
(2007) Modeling systems-level regulation of host immune responses
PLOS Computational Biology 3:e109.
https://doi.org/10.1371/journal.pcbi.0030109
- PubMed
- Google Scholar
(2009) Constraint-based network model of pathogen-immune system interactions
Journal of the Royal Society, Interface 6:599–612.
https://doi.org/10.1098/rsif.2008.0363
- PubMed
- Google Scholar
1. Thakar J
2. Pathak AK
3. Murphy L
4. Albert R
5. Cattadori IM
(2012) Network model of immune responses reveals key effectors to single and co-infection dynamics by a respiratory bacterium and a gastrointestinal helminth
PLOS Computational Biology 8:e1002345.
https://doi.org/10.1371/journal.pcbi.1002345
- PubMed
- Google Scholar
1. Vanalli C
2. Mari L
3. Righetto L
4. Casagrandi R
5. Gatto M
6. Cattadori IM
(2020) Within-host mechanisms of immune regulation explain the contrasting dynamics of two helminth species in both single and dual infections
PLOS Computational Biology 16:e1008438.
https://doi.org/10.1371/journal.pcbi.1008438
- PubMed
- Google Scholar
(2007) Comparative role of immunoglobulin A in protective immunity against the Bordetellae
Infection and Immunity 75:4416–4422.
https://doi.org/10.1128/IAI.00412-07
- PubMed
- Google Scholar
1. Woolfrey BF
2. Moody JA
(1991) Human infections associated with Bordetella bronchiseptica
Clinical Microbiology Reviews 4:243–255.
https://doi.org/10.1128/CMR.4.3.243
- PubMed
- Google Scholar

Decision letter

Niel Hens

Reviewing Editor; Hasselt University, Belgium
Aleksandra M Walczak

Senior Editor; CNRS LPENS, France
Samuel Alizon

Reviewer; Centre National de la Recherche Scientifique, France

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Gastrointestinal helminths increase Bordetella bronchiseptica shedding and host variation in supershedding" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Samuel Alizon (Reviewer #1).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) The authors present their results in the context of "supershedding" and, more generally, the idea that the variance in the number of secondary infections caused by an individual (i.e. R0) can matter in addition to the mean. Indeed, as popularized by Lloyd-Smith et al. (2005), even if the mean R0 is constant, an increased heterogeneity between individuals will affect disease emergence and spread. The authors appear to be presenting their work in this context and I agree that co-infections could be increasing the level of individual variation. However, in their results (and even more generally in the context of co-infections), I am not sure this is appropriate because it appears to me that the main consequence of the co-infection is to increase the mean of the distribution (rabbits shed more bacterial) rather than its heterogeneity (the proportion of rabbits that do not shed bacteria remains comparable with co-infections). If the authors really wish to keep the focus on the importance of heterogeneity, they should show that it matters and, for instance, estimate the heterogeneity parameter k of a negative binomial distribution (or a similar distribution) for the distribution of shedding rates (in Figure 3), assuming that this reflects individual R0. The other option is to focus on the mean value instead of heterogeneity, which I think can be done without loss for the integrity of the manuscript.

2) From Table 1, it appears that the authors are estimating at least 27 parameters using the model. This seems like a lot and the large confidence intervals make me wonder whether the model is identifiable. I am unsure this can be shown using likelihood profiles given the number of parameters but perhaps this could be studied by simulating datasets and calculating the mean relative error associated with the inference.

3) The authors mention model comparison in passing but given the number of parameters I think it might be worth exploring this in more details, especially since co-infections seem to be leading to similar patterns.

4) More is needed in the paragraph starting on line 309. The fact that neutrophils are produced during infections by both helminth species makes the explanation in lines 315-317 seem unconvincing. Why do the neutrophils have a different impact in the two helminth infections, and could it be related to different dynamics of helminth growth in the two species? The authors mention work in mice (Rolin et al.) that seems to show a different pattern of neutrophil dynamics and impact-should readers interpret that as merely differences between the immune responses of mice and rabbits, and if so, what differences are most likely? Possibly related to these points, it might be good to emphasize here that the model used did not account for differences in helminth infection intensity-would accounting for those differences in a future model be likely to shed light on the role of neutrophils?

5) More explanation would be helpful in a couple of places, especially regarding the dynamics in B. bronchiseptica singly infected hosts (lines 208-211). The text makes it sound like the data are just not informative as to where the peak is, but it seems clear that the peak occurred at or before the start of sampling. I was surprised that the peak was earlier for the single infections, where bacteria are thought to be replicating more slowly. Naively, I thought that with a slower replication rate, the peak would be later. How do the authors interpret that finding? Is it the case that the single infections are growing slower and also brought under control faster, resulting in an earlier peak than the double and triple infections.

6) The authors should pay particular attention to the specific comments raised as part of the public review and of course all remaining comments in a point-to-point reply.

Recommendations for the authors:

Reviewer #1:

lines 2-3: This sentence is a bit misleading because in Lloyd-Smith et al. the R0 is attributed to the individual.

Figure 1: The numbers on the top of each panel are unclear (I guess they refer to weeks?).

line 147-148: Is it an increase in shedding or in the likelihood to shed?

line 512-514: The writing is clear but I think the manuscript could gain in clarity by discussing a bit more within-host kinetics modelling, which is barely mentioned. Furthermore, regarding the model itself, modelling the immune response is more common now than it was 20 years ago, there are still many different models. For instance, here immune activation is assumed to depend only on parasite load. Referring to earlier models that made similar assumptions would help.

line 561-562: Do you need to make some assumptions regarding the independence of the variables to obtain the final likelihood function in equation 12?

Figure 6: It would be nice to also show in Figure 6 the within-host time series resulting from the model parameter inference and not only the "growth rates" (which are actually a bit unclear). This is already in Figure S1 but I think it would really improve the study.

Instead, showing likelihood profiles in the Appendix might be a good idea.

Reviewer #2:

Some modifications to the language (especially removing jargon and adding further explanation in places) would help make study accessible to a broader audience. For example, the phrase "rapid variation in individual shedding" (abstract) requires more explanation-perhaps "rapid temporal changes in individual shedding" would be clearer? Another example is line 96, where it seems strange to say that "events were null", and it might be more accessible to say that rabbits did not shed at many time points even though they were known to be infected and interacted with the petri dish (assuming I understood that correctly).

Lines 241-245 are confusing, starting with the phrase "negative bacterial growth". Assuming I'm interpreting the words in line 244 as a sharp decline in bacterial abundance, what does it mean to say that "the zero time to reach this peak were represented by the initial inoculum"?

I didn't get a lot out of Figures 7 and 8. The authors might consider moving these figures to the supplement, but regardless it would be helpful to include both parameter symbols and names/brief definitions on the x-axis labels or at the very least in the caption.

Reviewer #3:

Abstract, first two sentences. These sentences are a bit awkward. The authors should consider recasting them.

Abstract, "Model simulations revealed…". This makes no sense as written: simulations by themselves cannot tell us anything about the real world. Please revise this sentence to more accurately characterize the relationship between the data, the conclusions, and the model simulations.

Abstract, "…the rapid variation in individual shedding…". This sentence is very unclear.

Author Summary, line 2. Consider replacing "underline" with "underlie".

Author Summary, "experiments of rabbits together with mathematical modeling". It reads as though the experiments involved rabbits doing math!

Author Summary: "at the host level, …". This sentence is unclear. The authors should revisit it, and perhaps reconsider whether it belongs in the Author Summary.

When "type 1" and "type 2" are first introduced, they should be explained. At least, it should be made clear that "type 1" refers to Th1, etc.

ll 36ff. The connection between B. bronchiseptica and pertussis is tenuous and essentially irrelevant in this context. More generally, there is neither need nor value in this tangent.

Paragraph beginning on l 45. There are a number of facts mentioned here that are not obviously related to the authors' argument. The authors should consider whether these facts belong here. If they decide that they do belong, they should explain how.

ll 56ff. The authors describe some modeling work as if it were evidence. On the face of it, this is absurd. I recommend that they consider whether these sentences are needed, or contribute, to their study. If they conclude that they do, they should revise these sentences to put the earlier modeling work in context.

l 65. I find it strange that the authors jump from the experimental design immediately into the modeling without first describing the data that they generated.

l 72. Simulations cannot, by themselves, explain anything. Moreover, it requires great imagination to interpret the authors' model as "mechanistic".

l 80. "…every week or multiple times a week." This is most vague; the authors should be more precise.

l 91. Why was a nonparametric (Wilcoxon) test used? Does the result change if a parametric (e.g. t) test is used? Are there reasons to avoid such parametric tests? If so, what are they?

ll 93-95. Are the reported differences in median (?) shedding rate among the arms statistically significant? The broad and overlapping confidence intervals suggest not.

ll 105-111. This discussion is confusing and unclear.

ll 276-279. This is speculation, which is not in itself a problem, but it should be labeled as such.

ll 440-443. The choice of sampling times seems arbitrary. Can the authors describe the rationale behind these choices a bit more carefully.

l 449. "exemplifies" → "mimics".

ll 452-455. Explain the logistical constraints and technical difficulties.

ll 468, 470. This is an unwarranted conclusion. The authors should more carefully describe the conclusions that can be drawn from the cited study. In particular, the modeling study rests on strong assumptions about the underlying immunology.

l 568. "Weakly normal prior". "Weakness" is a relative term: the authors should describe the precise form of the priors they assume.

Figures:

Figure 1. There are several problems with this data visualization. First, it makes the "Alabama First" error, whereby the data are arranged according to an irrelevant variable (in this case, animal number within date-of-sacrifice). Second, the boxplots are not appropriate in many cases, since the data are too few or too non-normal. Third, the arbitrariness inherent in the log(1+CFU/s) metric makes it hard to interpret. I cannot confidently recommend any single visualization that will correct these problems: it will probably be necessary for the authors to experiment with, for example, simpler scatterplots, violin plots, and other approaches, before they find a more satisfactory plot or set of plots.

It may be that this figure is attempting to do too much. It seeks to convey information about the stereotypical time-course of infection and about the intra- and inter-animal variability, as well as the variation in both of these with coinfection. It might be helpful to design several figures that tackle each of the above individually.

Figure 2. This figure suffers from some of the same problems as Figure 1. In addition, there is too much cramping and overplotting to distinguish the individual-animal traces. The use of log(1+CFU/s) is problematic. Since the authors are using a zero-inflated model, it seems that the zeros don't belong here.

Finally, and at least as worryingly, the trends (smoothed curves) do not appear to represent the data at all. That is, the trends are not typical. Inasmuch as these median values are the point of contact with the models, this is quite problematic: it seems likely that even if the best-fitting models explain these averages well, they fail to represent any of the individual animals well.

Figure 3. Again, the log(1+CFU/s) metric is problematic. With a different, equally arbitrary, choice of time unit, the shape of these histograms might change appreciably.

ll 64ff. The authors appear to have made choices about the inclusion and exclusion of animals on an ad hoc basis. On its face, this raises questions about the reproducibility and reliability of their conclusions. However, it is strange that they have done so, since their zero-inflated model affords them a principled way of including all animals.

Figure S1 ought to be moved into the main text.

Figure 6. The model does not appear to do a good job in capturing the data. In particular, the individual traces and the overall trend appear to be biased downwards. This may be due to the presence of zeros in the data. If this is the case, then it is puzzling why the authors do not employ their zero-inflated model to focus attention on the non-zero data. Their choice instead to use log(1+CFU/s) is problematic in its own right as well, as discussed above.

This figure also suffers from the same problems as Figure 2: the individual data are not resolved and the individual model trajectories are too crowded.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Gastrointestinal helminths increase Bordetella bronchiseptica shedding and host variation in supershedding" for further consideration by eLife. Your revised article has been evaluated by Aleksandra Walczak (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

First set of remaining issues:

Figure 1-Could the authors include an x-axis label? I assume that would be something along the lines of rabbit ID.

Figures 2, 5, and 7: could the x-axis labels be fixed so that they are all legible? Also in Figure 5, the dark blue line is difficult to see against the thinner black lines.

The paragraph in lines 244-272 could use further proofreading since some phrases are difficult to parse, including "somehow represented by the infection dose" and "shedding was not statistically significant among the three groups".

Line 322-omit "prompted" or rephrase (meaning unclear).

Second set of remaining issues:

162 onwards and 335 – the definition of supershedding raises an interesting point – as noted here, supershedding is usually thought of as the maximum shedding however duration of shedding is an important component as well. Does it change the outcome at all, if the threshold was defined in terms of the integrated shedding profile?

172 Figure 2. I realise these are complex data to represent compactly, but I find it hard to follow the individual trajectories (e.g. to determine if single infections are consistently different from each other). Perhaps a few example trajectories in different colours might help?

Line 191 – Figure 3. It is hard to tell, but it looks like the negative binomial is a reasonable representation of the BT infections but that there may be greater systematic biases for the BG and BTG cases. Is this true and if so, statistically significant?

448 – I would be slightly more reserved about making a definitive statement on the impact of measured super-shedding in an experimental setting, an undefined use of 'contact' and the probability of it resulting in infection and/or a greater number of infections. There is a good chance they are right, but that's not quite the same as proving it.

684 and 932 onwards. Very good to see the convergence plots – it would be helpful also have the scale reduction factors, as its difficult see in the plots what the individual chains are actually doing – based on the text it should be fine, but having the results recorded would be useful.

902 – Good to see figures showing posterior estimates here – I much prefer them to tables. I do generally find it helpful to have figures showing the posterior distribution, preferably in correlation plots to show at least some of the interdependencies between parameters which is particularly important given the number of parameters involved. I realise it would add considerably to the appendix to show more detail for each individual, however, so long as it doesn't misrepresent the variation between individuals, I think an averaged correlation plot across all posteriors would be helpful.

https://doi.org/10.7554/eLife.70347.sa1

Author response

Essential revisions:

1) The authors present their results in the context of "supershedding" and, more generally, the idea that the variance in the number of secondary infections caused by an individual (i.e. R0) can matter in addition to the mean. Indeed, as popularized by Lloyd-Smith et al. (2005), even if the mean R0 is constant, an increased heterogeneity between individuals will affect disease emergence and spread. The authors appear to be presenting their work in this context and I agree that co-infections could be increasing the level of individual variation. However, in their results (and even more generally in the context of co-infections), I am not sure this is appropriate because it appears to me that the main consequence of the co-infection is to increase the mean of the distribution (rabbits shed more bacterial) rather than its heterogeneity (the proportion of rabbits that do not shed bacteria remains comparable with co-infections). If the authors really wish to keep the focus on the importance of heterogeneity, they should show that it matters and, for instance, estimate the heterogeneity parameter k of a negative binomial distribution (or a similar distribution) for the distribution of shedding rates (in Figure 3), assuming that this reflects individual R0. The other option is to focus on the mean value instead of heterogeneity, which I think can be done without loss for the integrity of the manuscript.

This is a very good point and important in the broader discussion on co-infection. We want to clarify that we do not quantify the proportion of hosts that shed or not shed, they all shed at some point during the experiment (the seven rabbits that never shed were excluded from the analysis).

This study investigated variation in the amount and frequency of shedding events at three levels: (i) between four types of infection, (ii) between hosts within each type of infection and iii- by every individual over time.

We have now clarified the objective of this work and in our updated results we show that helminths increase the level of bacteria shed, the frequency of shedding events (i.e. events when rabbits do not shed is 56% in the B group and 40%-25% in the co-infected groups) and variation in the magnitude of these events (Negative binomial k, BG: 0.37, BT: 0.45 and BTG: 0.21, no fitting for the B group because only three data-bins are available). This indicates that helminths affect both the mean of bacteria shed and variation in this mean, including the emergence of supershedding events.

We have revised figure 3 (the original plot was misleading) and we have now: i- used the same interval size (CFU/s=0.5) to classify shedding for the four infections, ii- provided the parameter k from fitting a Negative Binomial, iii- updated the 95th and 99th percentile thresholds as estimated using the whole data set (four infections together), and iv- presented the data as frequency of CFU/s events (i.e. not log-transformed as originally showed).

2) From Table 1, it appears that the authors are estimating at least 27 parameters using the model. This seems like a lot and the large confidence intervals make me wonder whether the model is identifiable. I am unsure this can be shown using likelihood profiles given the number of parameters but perhaps this could be studied by simulating datasets and calculating the mean relative error associated with the inference.

We are not sure how the Reviewer has calculated this number. The number of shared parameters within each infection group is 18 as given by equation (9) in the main text, while the number of individual parameters, given by equation (10), varies depending on the number of rabbits for each infection. To provide robust support to our approach, we have followed the advice and have now created new sections in Material and Methods (‘Model validation’) and Appendix (Appendix-4) where we investigate the ability of the model to recover the correct parameter values and show good model performance and accuracy for every estimated parameter.

3) The authors mention model comparison in passing but given the number of parameters I think it might be worth exploring this in more details, especially since co-infections seem to be leading to similar patterns.

We have now created a new section in Material and Methods (‘Model Selection’) and in the

Appendix (Appendix-5) where we compare three model formulations: the full model

(neutrophils+IgA+IgG) and two simplified versions (neutrophils+IgA and neutrophils+IgG). The full model was selected as the best compromise between goodness of fit and parsimony, and subsequently used in our analysis.

4) More is needed in the paragraph starting on line 309. The fact that neutrophils are produced during infections by both helminth species makes the explanation in lines 315-317 seem unconvincing. Why do the neutrophils have a different impact in the two helminth infections, and could it be related to different dynamics of helminth growth in the two species? The authors mention work in mice (Rolin et al.) that seems to show a different pattern of neutrophil dynamics and impact-should readers interpret that as merely differences between the immune responses of mice and rabbits, and if so, what differences are most likely? Possibly related to these points, it might be good to emphasize here that the model used did not account for differences in helminth infection intensity-would accounting for those differences in a future model be likely to shed light on the role of neutrophils?

We have carefully revised this section and removed some of the inconsistencies that caused confusion. We now discuss more coherently our neutrophil results and how our findings relate to previous studies in mice using Bordetella and other pathogens. To stress that we did not explicitly include the dynamics of the two helminths, we have included comments in the Introduction, Material and Methods and Results.

5) More explanation would be helpful in a couple of places, especially regarding the dynamics in B. bronchiseptica singly infected hosts (lines 208-211). The text makes it sound like the data are just not informative as to where the peak is, but it seems clear that the peak occurred at or before the start of sampling. I was surprised that the peak was earlier for the single infections, where bacteria are thought to be replicating more slowly. Naively, I thought that with a slower replication rate, the peak would be later. How do the authors interpret that finding? Is it the case that the single infections are growing slower and also brought under control faster, resulting in an earlier peak than the double and triple infections.

We have revised this section and provided more clarity on the dynamics of shedding in B. bronchiseptica only rabbits. Briefly, we note that model simulations place the peak right at the start of the trial, however, this should be taken with caution, since there are no shedding data to train the model during these first ten days (as noted in Material and Methods we used a different substrate for bacteria shedding that did not work). We also note that while we still report the model prediction at this early time this needs to be eventually validated in the laboratory.

Recommendations for the authors:

Reviewer #1:

lines 2-3: This sentence is a bit misleading because in Lloyd-Smith et al. the R0 is attributed to the individual.

The reviewer is correct in that Ro for microparasites is the number of secondary infections generated by a primary case; however, while this is the unit of measurement, essentially, Ro is a measurement of pathogen/parasite fitness and its ability to be transmitted.

Figure 1: The numbers on the top of each panel are unclear (I guess they refer to weeks?).

As originally noted in the legend: “The total number of rabbits monitored in each infection is reported in parenthesis at the top left of each plot”. To avoid confusion this detail has now been removed and the legend has been revised.

line 147-148: Is it an increase in shedding or in the likelihood to shed?

This is an increase in the probability of becoming supershedders (whether this is at the 95% or 99% percentile threshold), the sentence has been clarified.

line 512-514: The writing is clear but I think the manuscript could gain in clarity by discussing a bit more within-host kinetics modelling, which is barely mentioned. Furthermore, regarding the model itself, modelling the immune response is more common now than it was 20 years ago, there are still many different models. For instance, here immune activation is assumed to depend only on parasite load. Referring to earlier models that made similar assumptions would help.

We now briefly discuss previous models that investigated the within-host interaction between Bordetella and the immune response, and how our model differs from those; we also provide additional references of models that used assumptions similar to our work.

line 561-562: Do you need to make some assumptions regarding the independence of the variables to obtain the final likelihood function in equation 12?

We have clarified that there is the assumption of independence between the variables in the likelihood calculations. However, we also note that the variables have an interactive relationship that are captured in the dynamical models. Therefore, the assumption of independence inside the likelihood function is a plausible one.

Figure 6: It would be nice to also show in Figure 6 the within-host time series resulting from the model parameter inference and not only the "growth rates" (which are actually a bit unclear). This is already in Figure S1 but I think it would really improve the study.

Instead, showing likelihood profiles in the Appendix might be a good idea.

We think the reviewer refers to the original figure 5, the (original) figure 6 described bacterial shedding from the experimental and simulation data and included both individual and median results. We have updated figure 5 (current figure 6) and included the individual time series for the neutralization rates of the co-infection groups. We have not included individual time series for the growth and decay rates because of the very narrow CIs, which limit the ability to disentangle each individual trend. Simulations for the B group were done only at the group level. As advised by reviewer #3, we have now moved the original figure S1 in the main text (currently figure 5) to show the accuracy of model inference to experimental data.

Reviewer #2:

Some modifications to the language (especially removing jargon and adding further explanation in places) would help make study accessible to a broader audience. For example, the phrase "rapid variation in individual shedding" (abstract) requires more explanation-perhaps "rapid temporal changes in individual shedding" would be clearer? Another example is line 96, where it seems strange to say that "events were null", and it might be more accessible to say that rabbits did not shed at many time points even though they were known to be infected and interacted with the petri dish (assuming I understood that correctly).

We have revised the whole manuscript to improve clarify and accessibility to a broader audience.

Lines 241-245 are confusing, starting with the phrase "negative bacterial growth". Assuming I'm interpreting the words in line 244 as a sharp decline in bacterial abundance, what does it mean to say that "the zero time to reach this peak were represented by the initial inoculum"?

We agree that this oxymoron is confusing and have revised the paragraph.

I didn't get a lot out of Figures 7 and 8. The authors might consider moving these figures to the supplement, but regardless it would be helpful to include both parameter symbols and names/brief definitions on the x-axis labels or at the very least in the caption.

We have followed the Reviewer’s suggestion and moved these figures in the Appendix (now Appendix 2- figures 1-2); we also revised the labeling as advised to increase clarity.

Reviewer #3:

Abstract, first two sentences. These sentences are a bit awkward. The authors should consider recasting them.

Abstract, "Model simulations revealed…". This makes no sense as written: simulations by themselves cannot tell us anything about the real world. Please revise this sentence to more accurately characterize the relationship between the data, the conclusions, and the model simulations.

Abstract, "…the rapid variation in individual shedding…". This sentence is very unclear.

The abstract has been revised to improve clarity.

Author Summary, line 2. Consider replacing "underline" with "underlie".

Author Summary, "experiments of rabbits together with mathematical modeling". It reads as though the experiments involved rabbits doing math!

Author Summary: "at the host level, …". This sentence is unclear. The authors should revisit it, and perhaps reconsider whether it belongs in the Author Summary.

Accordingly and as noted above, the summary has been revised.

When "type 1" and "type 2" are first introduced, they should be explained. At least, it should be made clear that "type 1" refers to Th1, etc.

We have replaced type 1 and type 2 with Th1 and Th2, which are terms more commonly known. When first mentioned, we refer to ‘Th1 inflammatory response’ or ‘Th2 anti-inflammatory response’ to provide more information on the type of response.

ll 36ff. The connection between B. bronchiseptica and pertussis is tenuous and essentially irrelevant in this context. More generally, there is neither need nor value in this tangent.

We included this connection as relevant in the context of understanding the dynamics of shedding of Bordetella species, including in humans where whooping cough remains an infection of serious concern. In this context, the connection with humans gives broad relevance to the study beyond the rabbit system.

Paragraph beginning on l 45. There are a number of facts mentioned here that are not obviously related to the authors' argument. The authors should consider whether these facts belong here. If they decide that they do belong, they should explain how.

This section discusses co-infection and provides a background on what has been previously done on co-infections with B. bronchiseptica and helminths. This is also important to understand the system and the rationale of our study, and to prepare the reader to the following section on the objective and general working approach. We have revised the Introduction and this section accordingly.

ll 56ff. The authors describe some modeling work as if it were evidence. On the face of it, this is absurd. I recommend that they consider whether these sentences are needed, or contribute, to their study. If they conclude that they do, they should revise these sentences to put the earlier modeling work in context.

In this section we describe findings from our previous work using experimental data and modeling on the same system that we use in the current study. This is relevant for two reasons: first, it provides useful information to understand the system and second, it informs on a possible mechanism of immune regulation of Bordetella in rabbits. This section has been revised and references on similar work in mice added. Please, see additional previous comment regarding the need of this whole paragraph.

l 65. I find it strange that the authors jump from the experimental design immediately into the modeling without first describing the data that they generated.

We have now revised this section to avoid the ‘jump from experiments into modeling’ and provided more information on the data used, including referring to previous work.

l 72. Simulations cannot, by themselves, explain anything. Moreover, it requires great imagination to interpret the authors' model as "mechanistic".

We have revised the sentence accordingly, however, we disagree with the reviewer regarding the term ‘mechanistic’ as used in this context. The model we propose offers a parsimonious, and simplified, explanation of the mechanism in which the three immune variables affect the temporal dynamics of Bordetella infection and shedding (please, see further comments below).

l 80. "…every week or multiple times a week." This is most vague; the authors should be more precise.

We have now clarified that shedding data were collected every week or two-three times a week for every rabbit and refer to Material and Methods (Bacteria shed enumeration) for additional details.

l 91. Why was a nonparametric (Wilcoxon) test used? Does the result change if a parametric (e.g. t) test is used? Are there reasons to avoid such parametric tests? If so, what are they?

Thank you for pointing this out. We found a make a mistake in the initial choice of the hypothesis test. The appropriate test to use is the Kruskal Wallis test, which we used to compare the group medians among the four types of infection with Bonferroni’s correction for multiple testing. We choose a non-parametric test because the data are highly aggregated (see Figure 3).

ll 93-95. Are the reported differences in median (?) shedding rate among the arms statistically significant? The broad and overlapping confidence intervals suggest not.

The reviewer is correct, we found no differences between the co-infected groups and included the significant term (p>0.05) of our analysis to stress our point.

ll 105-111. This discussion is confusing and unclear.

This section has been revised to avoid confusion.

ll 276-279. This is speculation, which is not in itself a problem, but it should be labeled as such.

We have revised the section to avoid definitive conclusions.

ll 440-443. The choice of sampling times seems arbitrary. Can the authors describe the rationale behind these choices a bit more carefully.

The experimental design is driven by the dynamics/life cycles of the three infectious agents. Specifically, we followed the standard design of Bordetella infection in mice, which considers important phases in the bacteria-host interaction (e.g. Kirimanjeswara et al. 2003). We also included additional time points to increase accuracy in our dataset and representative phases in the life cycle of the two helminths (Murphy et al. 2022, 2013, Cattadori et al. 2019) that might be important in affecting Bordetella. We have revised the sentence and included appropriate references to support our sampling design.

l 449. "exemplifies" → "mimics".

We made the change as suggested.

ll 452-455. Explain the logistical constraints and technical difficulties.

We consider this information irrelevant to the Methods, however, we have positively addressed this comment as follow: “The use of a different number of hosts and frequency of sampling was determined by logistical constraints (i.e. personnel availability).”

ll 468, 470. This is an unwarranted conclusion. The authors should more carefully describe the conclusions that can be drawn from the cited study. In particular, the modeling study rests on strong assumptions about the underlying immunology.

While we do agree that ‘the modeling study rests on strong assumptions about the underlying immunology’, the cited models do indeed build on data from laboratory experiments and current knowledge on the host immune response to Bordetella infection. Our modeling assumptions follow a similar rationale using data from our laboratory experiments. We have revised the sentence to avoid definitive statements.

l 568. "Weakly normal prior". "Weakness" is a relative term: the authors should describe the precise form of the priors they assume.

We have now provided the prior values in table 2 and noted that they were generated from a Normal distribution.

Figures:

Figure 1. There are several problems with this data visualization. First, it makes the "Alabama First" error, whereby the data are arranged according to an irrelevant variable (in this case, animal number within date-of-sacrifice). Second, the boxplots are not appropriate in many cases, since the data are too few or too non-normal. Third, the arbitrariness inherent in the log(1+CFU/s) metric makes it hard to interpret. I cannot confidently recommend any single visualization that will correct these problems: it will probably be necessary for the authors to experiment with, for example, simpler scatterplots, violin plots, and other approaches, before they find a more satisfactory plot or set of plots.

It may be that this figure is attempting to do too much. It seeks to convey information about the stereotypical time-course of infection and about the intra- and inter-animal variability, as well as the variation in both of these with coinfection. It might be helpful to design several figures that tackle each of the above individually.

Figure 1 has been designed to mimic our experiential design, a stereotypical visualization yet an effective way to convey a simple message while avoiding graphic redundancy. Within each infection rabbits are ordered from left to right following the experimental time. We sacrificed groups of 4 rabbits at a time, groups sampled early on the left and groups sampled late on the right; to facilitate the vision of the time course we have reported the week of sampling for each group of animals. The order of these four rabbits within each group is arbitrary since this is the only time we report the ID of the host. We tried other approaches and the boxplot was still the best compromise between a representation of the experimental design, which we wanted to provide, and general information on the type of infections and metrics of each host. As previously reported, the log(CFU/s+1) is a useful choice for the visualization of highly variable data.

Figure 2. This figure suffers from some of the same problems as Figure 1. In addition, there is too much cramping and overplotting to distinguish the individual-animal traces. The use of log(1+CFU/s) is problematic. Since the authors are using a zero-inflated model, it seems that the zeros don't belong here.

Finally, and at least as worryingly, the trends (smoothed curves) do not appear to represent the data at all. That is, the trends are not typical. Inasmuch as these median values are the point of contact with the models, this is quite problematic: it seems likely that even if the best-fitting models explain these averages well, they fail to represent any of the individual animals well.

Figure 2 represents the experimental longitudinal data on shedding and highlights the high heterogeneity in the shedding events both within and between hosts. The smoothed individual curves are challenged by the high variation between consecutive time points, this is also highlighted in the median trend. To help with the visualization, data were transformed to log(CFU/s+1), as mentioned above. We have revised this figure and now plot the individual trajectories by simply joining the events with straight segments but keeping the smoothed median fitted to the group data. We have also revised the legend and our comments to this figure in the main text. We tried alternative plotting but visualization did not improve, given the nature of the data and what we wanted to convey.

Figure 3. Again, the log(1+CFU/s) metric is problematic. With a different, equally arbitrary, choice of time unit, the shape of these histograms might change appreciably.

The frequency distribution is now based on CFU/s data (not log-transformed) with intensity of shedding grouped in classes of 0.5 unit intervals (these are not time intervals), plus the initial class 0 for the ‘null-shedding’. This approach has been applied to each type of infection. We have also calculated the 99th and 95th percentile thresholds from the whole data set and used these cut-offs to estimate the number of supershedding events/hosts in each infection group. Finally, we have fitted a Negative Binomial to the frequency distribution. This updated figure provides new results that are reported and discussed in the main text.

ll 64ff. The authors appear to have made choices about the inclusion and exclusion of animals on an ad hoc basis. On its face, this raises questions about the reproducibility and reliability of their conclusions. However, it is strange that they have done so, since their zero-inflated model affords them a principled way of including all animals.

The reviewer probably meant l87 of the original manuscript. To avoid confusion, we have now clarified that since we were interested in the dynamics of shedding, rabbits that never shed (n = 7) were excluded from the subsequent analysis. This number of rabbits is also too small to perform any meaningful analysis on the immune response associated with rabbits that do not shed. Please see detailed comments on the use of animals above.

Figure S1 ought to be moved into the main text.

The figure has now been moved into the main text as suggested.

Figure 6. The model does not appear to do a good job in capturing the data. In particular, the individual traces and the overall trend appear to be biased downwards. This may be due to the presence of zeros in the data. If this is the case, then it is puzzling why the authors do not employ their zero-inflated model to focus attention on the non-zero data. Their choice instead to use log(1+CFU/s) is problematic in its own right as well, as discussed above.

This figure also suffers from the same problems as Figure 2: the individual data are not resolved and the individual model trajectories are too crowded.

We have revised the quality of figure 6 (now figure 7) to make our results clearer. We also realized that our comments were not completely accurate and, thus, the main text has now been revised. Briefly, we note that “simulations of individual time series captured relatively well the general trend and the high level of shedding early in the infection, particularly for BTG and BT, and less clearly for BG”. By changing the individual trajectories to a black tone, we can now see the general good fit to the data.

We did apply the zero-inflated relationship between shedding and infection and used the logtransformed data to reduce variation and skewed trend. Specifically, shedding data were modeled with a zero-inflated distribution where the 0 values constitute the zero-inflated part and the non-zero values are captured by a log-Normal distribution. Please, see our previous comments on the use of log-transformed data.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

First set of remaining issues:

Figure 1-Could the authors include an x-axis label? I assume that would be something along the lines of rabbit ID.

Good point, we have now included the x-label “rabbit ID’ as suggested.

Figures 2, 5, and 7: could the x-axis labels be fixed so that they are all legible? Also in Figure 5, the dark blue line is difficult to see against the thinner black lines.

We have revised the figures and provided a clearer format by adjusting the x-axis and making the contrast between lines clearer.

The paragraph in lines 244-272 could use further proofreading since some phrases are difficult to parse, including "somehow represented by the infection dose" and "shedding was not statistically significant among the three groups".

We have revised the section to avoid confusion.

Line 322-omit "prompted" or rephrase (meaning unclear).

The sentence has been revised accordingly.

Second set of remaining issues:

162 onwards and 335 – the definition of supershedding raises an interesting point – as noted here, supershedding is usually thought of as the maximum shedding however duration of shedding is an important component as well. Does it change the outcome at all, if the threshold was defined in terms of the integrated shedding profile?

This is an interesting consideration. We have further investigated the shedding data by calculating the area under the curve for each individual trajectory as composite trapezoids of the area under the linear lines connecting all the shedding values and the x-axis. We found that at a 95% threshold (defined from all infections) the supershedding events were as follow, B: 0, BG: 1 (4.5%), BT: 0 and BTG: 3 (12.5%), while at a 99% only the BG group (1, 4.5%) continued to have supershedding events. The high BG supershedders seem to substantially drive this pattern and affect the other groups. While interesting, we did not include this result in the manuscript to avoid confusion with the more classical approach we present and the general message we wanted to convey.

172 Figure 2. I realise these are complex data to represent compactly, but I find it hard to follow the individual trajectories (e.g. to determine if single infections are consistently different from each other). Perhaps a few example trajectories in different colours might help?

We have revised the figure and made it in a bigger size, we have also highlighted the trajectories of a few individuals for every infection, as suggested.

Line 191 – Figure 3. It is hard to tell, but it looks like the negative binomial is a reasonable representation of the BT infections but that there may be greater systematic biases for the BG and BTG cases. Is this true and if so, statistically significant?

We have now noted that the negative binomial distribution was not significantly different from the frequency of each co-infection dataset (p>0.05). In the original version we also provided the CIs of the aggregation parameter k.

448 – I would be slightly more reserved about making a definitive statement on the impact of measured super-shedding in an experimental setting, an undefined use of 'contact' and the probability of it resulting in infection and/or a greater number of infections. There is a good chance they are right, but that's not quite the same as proving it.

We have revised the section to avoid an overstatement in our reasoning.

684 and 932 onwards. Very good to see the convergence plots – it would be helpful also have the scale reduction factors, as its difficult see in the plots what the individual chains are actually doing – based on the text it should be fine, but having the results recorded would be useful.

We have now included Table 1 in Appendix 3 where we report the scale reduction factors of each parameter for the four types of infections.

902 – Good to see figures showing posterior estimates here – I much prefer them to tables. I do generally find it helpful to have figures showing the posterior distribution, preferably in correlation plots to show at least some of the interdependencies between parameters which is particularly important given the number of parameters involved. I realise it would add considerably to the appendix to show more detail for each individual, however, so long as it doesn't misrepresent the variation between individuals, I think an averaged correlation plot across all posteriors would be helpful.

We have included in Appendix 3 the correlations and related plots of the posteriors for every infection group, as advised.

https://doi.org/10.7554/eLife.70347.sa2

Article and author information

Author details

Nhat TD Nguyen
1. Center for Infectious Disease Dynamics, The Pennsylvania State University, University Park, United States
2. Department of Biology, The Pennsylvania State University, University Park, United States
Contribution
Formal analysis, Validation, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Ashutosh K Pathak

Competing interests
No competing interests declared
Ashutosh K Pathak
1. Center for Infectious Disease Dynamics, The Pennsylvania State University, University Park, United States
2. Department of Biology, The Pennsylvania State University, University Park, United States
3. Department of Infectious Diseases, University of Georgia, Athens, United States
Contribution
Data curation, Methodology, Writing – review and editing, Animal work, laboratory analysis

Contributed equally with
Nhat TD Nguyen

Competing interests
No competing interests declared
Isabella M Cattadori
1. Center for Infectious Disease Dynamics, The Pennsylvania State University, University Park, United States
2. Department of Biology, The Pennsylvania State University, University Park, United States
Contribution
Conceptualization, Funding acquisition, Investigation, Writing – original draft, Project administration, Writing – review and editing

For correspondence
imc3@psu.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6618-316X

Funding

Human Frontier Science Program (RGP0020/2007-C)

Ashutosh K Pathak

National Science Foundation (1145697)

Ashutosh K Pathak

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors thank Kathleen Creppage and Chad Pelensky for their research assistance during animal work and sample processing. This study, IMC and AKP were supported by Human Frontier Science Program (RGP0020/2007 C) and National Science Foundation-DEB (1145697). The funders had no role in study design, data collection and analysis, manuscript preparation and decision to publish.

Ethics

Animals were housed in individual cages with food and water ad libitum and a 12hr day/night cycle, in compliance with Animal Welfare Act regulations as well as the Guide for the Care and Use of Laboratory Animals. All animal procedures, including infections with the bacterium and the two helminth species, weekly blood collection and pathogen/parasite sampling at fixed time, were approved by the Institutional Animal Care and Use Committee of The Pennsylvania State University (IACUC 26082). All animal work complied with guidelines as reported in the Guide for the Care and Use of Laboratory Animals. 8th ed. National Research Council of the National Academies, National Academies Press Washington DC.

Senior Editor

Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

Niel Hens, Hasselt University, Belgium

Reviewer

Samuel Alizon, Centre National de la Recherche Scientifique, France

Publication history

Preprint posted: May 6, 2021 (view preprint)
Received: May 14, 2021
Accepted: October 11, 2022
Version of Record published: November 8, 2022 (version 1)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

218

Page views
30

Downloads
0

Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Nhat TD Nguyen
Ashutosh K Pathak
Isabella M Cattadori

(2022)

Gastrointestinal helminths increase Bordetella bronchiseptica shedding and host variation in supershedding

eLife 11:e70347.

https://doi.org/10.7554/eLife.70347

Share this article

Cite this article

Observed levels of B. bronchiseptica shedding by contact of infected rabbits with a Bordet-Gengou agar petri dish.

Observed level of B.bronchiseptica shedding by days post infection (Dpi) and type of infection.

Frequency of observed B.bronchiseptica shedding events for the four types of infection.

B. bronchiseptica shedding on BG-blood agar petri dishes.

Estimated immune responses in the blood by time and type of infection.

Estimated rates of infection in the whole respiratory tract and immune response in the blood by time and type of infection.

Description and unit of parameters and variables used in the dynamical and observed models.

Model prior and posterior parameters for the immune response and B. bronchiseptica infection.

Estimated dynamics of B. bronchiseptica shedding by time and type of infection.

Parameter estimates from posterior distribution for BG.

Parameter estimates from posterior distribution for BT.

Parameter estimates from posterior distribution for BTG.

Relationship between B.bronchiseptica peak of shedding and percentile changes in neutrophils, IgA, and IgG.

Relationship between B.bronchiseptica time to peak of shedding and percentile changes in neutrophils, IgA, and IgG.

MCMC trace plot for B.

MCMC trace plot for BG.

MCMC trace plot for BT.

MCMC trace plot for BTG.

MCMC posteriors from B.

MCMC posteriors from BG.

MCMC posteriors from BT.

MCMC posteriors from BTG.

Scale reduction factors for each parameter of the four types of infection.

Model validation using group-level parameters.

Bayesian Information Criterion (BIC) and level of model complexity (N), including pair-wise Δ BIC with the best fitted model.

Author details

Nhat TD Nguyen

Contribution

Contributed equally with

Competing interests

Ashutosh K Pathak

Contribution

Contributed equally with

Competing interests

Isabella M Cattadori

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading